acer

29759 Reputation

29 Badges

19 years, 316 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

The LinearSolve command in the LinearAlgebra package can handle this.

Let's write out the coefficients (of x) for your polynomials, as column Vectors. Form a Matrix of the ones from S.

Then call LinearSolve.

restart;

with(LinearAlgebra,LinearSolve):

 

S:={x+4, 3*x-7}: f:=-5*x+10:

T := max(degree~(S,x))+1:

A:=Matrix(T,nops(S),(i,j)->coeff(S[j],x,i-1));

Matrix(2, 2, {(1, 1) = 4, (1, 2) = -7, (2, 1) = 1, (2, 2) = 3})

B:=Vector(T,(i)->coeff(f,x,i-1));

Vector(2, {(1) = 10, (2) = -5})

LinearSolve(A,B);

Vector(2, {(1) = -5/19, (2) = -30/19})

 

SC := {x+4, 2*x+8}:

T := max(degree~(SC,x))+1:

AC:=Matrix(T,nops(SC),(i,j)->coeff(SC[j],x,i-1));

Matrix(2, 2, {(1, 1) = 4, (1, 2) = 8, (2, 1) = 1, (2, 2) = 2})

LinearSolve(AC,B);

Error, (in LinearAlgebra:-LinearSolve) inconsistent system

 

S:={x^2+x+4, 7*x^2+3*x, x+1}: f:=-5*x+10:

T := max(degree~(S,x))+1:

A:=Matrix(T,nops(S),(i,j)->coeff(S[j],x,i-1));

Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 4, (2, 1) = 1, (2, 2) = 3, (2, 3) = 1, (3, 1) = 0, (3, 2) = 7, (3, 3) = 1})

B:=Vector(T,(i)->coeff(f,x,i-1));

Vector(3, {(1) = 10, (2) = -5, (3) = 0})

LinearSolve(A,B);

Vector(3, {(1) = -15/2, (2) = -5/8, (3) = 35/8})

Download LS_ex2.mw

note: You wrote, "corresponding linear system", which indicates that you recognize that this problem can be expressed as a linear system. Using LinearSolve might then not be unexpected. Reformulating the polynomial relationships as a linear system, programmatically, might then be useful for you.

Notice that you were using the %d format code for this item, in printf. That's for integers. But you mentioned wanting them as 1.0,1.2,...2.0, which are floats. So you could use the %2.1f format, in printf.

Now, you already have the float values you mentioned wanting, available directly via t[i] .

If you didn't already have them using t[i] then you could also get n values evenly spaced between a and b by using the formula,
   a + (b-a)*(i-1)/(n-1)
which for this example would be,
   1 + (2-1)*(i-1)/(n-1)
which can be shortened to,
   1 + (i-1)/(n-1)

restart;

actual:=<2,3.5136,5.8016,9.1136,13.7376,20>:
t:=<1,1.2,1.4,1.6,1.8,2>:
n:=6:

P:=Statistics:-LinearFit([1,x],t,actual,x):

printf("t    Actual P  Model P  Error\n");
for i from 1 to n do

   model:=subs(x=t[i],P):
   err:=actual[i]-model:

   #printf("%2.1f %7.4f% 10.4f% 9.4f\n",
   #        1+(i-1)/(n-1),actual[i],model,err);

   printf("%2.1f %7.4f% 10.4f% 9.4f\n",
          t[i],actual[i],model,err);

end do:

t    Actual P  Model P  Error
1.0  2.0000    0.1717   1.8283
1.2  3.5136    3.7141  -0.2005
1.4  5.8016    7.2565  -1.4549
1.6  9.1136   10.7989  -1.6853
1.8 13.7376   14.3413  -0.6037
2.0 20.0000   17.8837   2.1163

Download frm.mw


For future programming, you might also care to think about how the following work. These kinds of construction occur pretty often in programming. It can be a useful thing to get familiar with.

<seq(1.0..2.0,0.2)>;
n:=6:
<seq(1.0..2.0,numelems=n)>;
<seq(1.0..2.0,(2-1)/(n-1))>;
<seq(1.0+(i-1)/(n-1),i=1..n)>;
Vector(n,i->1.0+(i-1)/(n-1));

You are using module locals like Chevron2A as the keys of some of your tables.

That will prevent you from accessing table entries if you subsequently try and use (say) the global names like Chevron2A as the key, trying that access outside the module.

You could, instead, use global name :-Chevron2A as the table key when creating the entry. Then you can use that as key, for access elsewhere. [edit] I would prefer this, over a module export. Using strings is also a good possibility, as Carl has mentioned after I posted this Answer. If you use the global name as key then you ought to quote it (and protect from evaluation), in case it is also assigned some value, eg. ':-Chevon2A'.

Or you could instead declare Chevron2A as an export of the module. Then, later, you could use it as key in its long-form (fully qualified name),  calculateblahcubitblah:-Chevron2A. Or you could make the module a package, so after doing with(calculateblahcubitblah) its exports become the new bindings for reference via the short-form Chevron2A. [edit] This is not my favorite choice, and I can only think of convoluted scenarios in which it'd provide extra contextual usefulness.

The code snippets that you showed above indicate that this problem that I've described occurs in your code. I cannot tell whether you also have other issues, eg. last-name-eval of tables assigned to names. I suggest using eval(nameoftable) when putting into lists or if returning them straight from procedures.

Some of your procedures will throw an error for certains ranges of inputs alpha and delta.

The plots:-inequal command appears to get slightly confused by that situation, and in consequence can misrepresent the feasible region..

Below and use short wrappers for those three procedures, so that when failing to produce a result (throwing an error) they instead return either a large positive or large negative value. Those can be used appropriately for the inequality tests, so that a more clear/robust representation of the feasible region is attained.

Double check, for correctness.

restartNULL

with(plots); with(RealDomain)

NULL

c := 1; cr := 0.3e-1*c; u := 1; sExp := 0.6e-1*c; s := .65*c; v := 3*cNULL

NULL

FirmModelPP := proc (alpha, beta, delta) local p0, xi0, q0, Firmpf0, G01, Recpf0, Unsold0, Environ0; option remember; xi0 := 1; p0 := min(s+sqrt((v-s)*(c-s)), delta*v+sExp); q0 := u*(v-p0)/(v-s); f(N) := 1/u; F(N) := N/u; G01 := int(F(N), N = 0 .. q0); Firmpf0 := (p0-c)*q0-(p0-s)*G01; Recpf0 := (sExp-cr)*xi0*q0; Unsold0 := G01; Environ0 := q0+Unsold0; return p0, q0, Firmpf0, Recpf0, Environ0, Unsold0 end proc

NULL

NULLNULLNULL

FirmModelFC := proc (alpha, beta, delta) local p00, xi00, q00, Firmpf00, G001, G002, Recpf00, Unsold00, Environ00, pr00; option remember; xi00 := 1; p00 := s+sqrt((v-s)*(c-s)); q00 := alpha*u*(v-p00)/(v-s); f(N) := 1/u; F(N) := N/u; G001 := int(F(N), N = 0 .. q00/alpha); G002 := int(F(N), N = 0 .. beta*xi00*q00/(1-alpha)); pr00 := p00-delta*v; Firmpf00 := (p00-c)*q00-alpha*(p00-s)*G001; Recpf00 := xi00*q00*(sExp-cr)+(pr00-sExp)*(beta*xi00*q00-(1-alpha)*G002); Unsold00 := alpha*G001; Environ00 := q00+Unsold00; return p00, q00, Firmpf00, Recpf00, Environ00, Unsold00 end proc

NULLNULLNULL

NULLNULL

NULL

FirmModelHmax := proc (alpha, beta, delta) local q, p, pr, FirmpfSiS, F1, G1, G2, G3, RecpfSiS, sol, UnsoldSiS, EnvironSiS, p0, OldSoldPrim, xi, h, ps, qs, prs, prof1m, prof2m; option remember; xi := 1; prs := ps-delta*v; prof1m := (ps-c)*qs+((1/2)*beta^2*xi^2*qs^2/(u*(1-alpha))-(1/2)*(1+beta*xi)^2*qs^2/u)*(ps-s)+(prs-sExp)*(beta*xi*qs-(1/2)*beta^2*xi^2*qs^2/(u*(1-alpha))); prof2m := (ps-c)*qs-(1/2)*(ps-s)*qs^2/(alpha*u)+(prs-sExp)*(beta*xi*qs-(1/2)*beta^2*xi^2*qs^2/(u*(1-alpha))); if alpha <= 1/(1+beta*xi) then p, q := (eval([ps, qs], solve({diff(prof1m, qs) = 0, qs = alpha*u*(v-ps)/(v-s), c < ps, sExp+delta*v < ps, qs < 2*u*(1-alpha)/(beta*xi)}, [ps, qs])[1]))[]; `h&Assign;`*(p-delta*v-sExp)/(p-delta*v); FirmpfSiS := eval(prof1m, [ps = p, qs = q, prs = p-delta*v]); RecpfSiS := (sExp-cr)*xi*q; UnsoldSiS := (1/2)*(1+beta*xi)^2*q^2/u-(1/2)*beta^2*xi^2*q^2/(u*(1-alpha)); EnvironSiS := q+UnsoldSiS else p, q := (eval([ps, qs], solve({diff(prof2m, qs) = 0, qs = alpha*u*(v-ps)/(v-s), c < ps, sExp+delta*v < ps, qs < 2*u*(1-alpha)/(beta*xi)}, [ps, qs])[1]))[]; h := (p-delta*v-sExp)/(p-delta*v); FirmpfSiS := eval(prof2m, [ps = p, qs = q, prs = p-delta*v]); RecpfSiS := (sExp-cr)*xi*q; UnsoldSiS := (1/2)*q^2/(alpha*u); EnvironSiS := q+UnsoldSiS end if; return p, q, FirmpfSiS, RecpfSiS, EnvironSiS, UnsoldSiS, h, OldSoldPrim, xi end proc

NULL

NULLNULL

pltPP3A := plot('FirmModelPP(alpha, .35, .40)[3]', alpha = 0. .. 1.0, color = red, legend = "", style = pointline, labels = [alpha, "Firm Profit"], labeldirections = ["horizontal", "vertical"], symbolsize = 10, axes = boxed, symbol = box, numpoints = 10, adaptive = false, thickness = 1.0, view = [0 .. 1, 0 .. .18])

NULL

pltPP3B := plot([[0., eval('FirmModelPP(0., .35, .40)[3]', alpha = 0.)]], color = red, legend = ["PP"], style = point, symbol = box, symbolsize = 10, axes = boxed, view = [0 .. 1, 0 .. .18])

NULL

pltFC3A := plot('FirmModelFC(alpha, .35, .40)[3]', alpha = 0. .. 1, linestyle = [solid], color = black, legend = "FC", labels = [alpha, "Firm Profit"], labeldirections = ["horizontal", "vertical"], axes = boxed, adaptive = false, thickness = .7)

NULL

NULLNULLNULLNULLNULL

pltHmax3A := plot('FirmModelHmax(alpha, .35, .40)[3]', alpha = 0. .. 1.0, linestyle = [dashdot], color = brown, legend = [SiS(h__max)], labels = [alpha, "Firm Profit"], labeldirections = ["horizontal", "vertical"], symbolsize = 10, numpoints = 50, adaptive = false, thickness = 1.0, axes = boxed)

NULL

#H:=plot3d((a,d)->FirmModelHmax(a, 0.35, d)[3]-FirmModelFC(a, 0.35, d)[3],
#       0..1,0..0.5,grid=[49,49],style=surface,orientation=[-90,0,0],lightmodel=none,
#       color=((a,d)->1/3*(1+signum(FirmModelHmax(a, 0.35, d)[3]-FirmModelFC(a, 0.35, d)[3]))),
#       adaptmesh=true
#      );

FMHmax := proc (a, d, lowhigh::(identical("low", "high"))) local res; option remember; if not [a, d]::(list(numeric)) then return 'procname(args)' end if; try res := FirmModelHmax(a, .35, d)[3] catch: res := undefined end try; if res::numeric then res elif lowhigh = "low" then -199 else 199 end if end proc; FMFC := proc (a, d, lowhigh::(identical("low", "high"))) local res; option remember; if not [a, d]::(list(numeric)) then return 'procname(args)' end if; try res := FirmModelFC(a, .35, d)[3] catch: res := undefined end try; if res::numeric then res elif lowhigh = "low" then -199 else 199 end if end proc; FMPP := proc (a, d, lowhigh::(identical("low", "high"))) local res; option remember; if not [a, d]::(list(numeric)) then return 'procname(args)' end if; try res := FirmModelPP(a, .35, d)[3] catch: res := undefined end try; if res::numeric then res elif lowhigh = "low" then -199 else 199 end if end proc

 

 

P1 := plots:-inequal(FMFC(a, d, "high") <= FMHmax(a, d, "low"), a = 0 .. 1, d = 0 .. .5, optionsfeasible = [[color = "Spring 1"]])
NULL

P2 := plots:-inequal(FMFC(a, d, "low") >= FMHmax(a, d, "high"), a = 0 .. 1, d = 0 .. .5, optionsfeasible = [[color = "Spring 3"]])

plots:-display(P1, P2)

P3 := plots:-inequal(FMPP(a, d, "high") <= FMHmax(a, d, "low"), a = 0 .. 1, d = 0 .. .5, optionsfeasible = [[color = "Spring 2"]])
NULL

P4 := plots:-inequal(FMPP(a, d, "low") >= FMHmax(a, d, "high"), a = 0 .. 1, d = 0 .. .5, optionsfeasible = [[color = "Spring 4"]])

plots:-display(P3, P4)

Download ConflcitInequal_acc.mw

These results for n=2..15 are computed exactly.

restart;

str:="1/(3 n Sqrt[Pi] Gamma[2 - n] Gamma[2 + n])*2^(-7 - 4 n) (3 256^
  n Gamma[2 - n] Gamma[-(1/2) +
    n] (8 n (-3 + 5 n + 2 n^2) Hypergeometric2F1[1 - n,
      2 - n, -2 n, 3/
      4] - (-1 +
       n) (4 (-4 + 7 n + 2 n^2) Hypergeometric2F1[2 - n, 2 - n,
         1 - 2 n, 3/4] - (-10 + 3 n + n^2) Hypergeometric2F1[
         2 - n, 3 - n, 2 - 2 n, 3/4])) -
 9^n Gamma[-(1/2) - n] Gamma[
   2 + n] (8 n (5 + 11 n + 2 n^2) Hypergeometric2F1[1 + n, 2 + n,
      2 n, 3/4] -
    3 (1 + n) (4 (4 + 9 n + 2 n^2) Hypergeometric2F1[2 + n, 2 + n,
          1 + 2 n, 3/4] -
       3 (6 + 5 n + n^2) Hypergeometric2F1[2 + n, 3 + n,
         2 (1 + n), 3/4])))":

ee := unapply(simplify(MmaTranslator:-FromMma(str)),n);

proc (n) options operator, arrow; -(1/24)*(-48*2^(-4+4*n)*(n+1)*GAMMA(3/2+n)^2*cos(Pi*n)*(n-1/2)*hypergeom([1-n, 1-n], [1-2*n], 3/4)-6*2^(-4+4*n)*GAMMA(3/2+n)^2*cos(Pi*n)*(n-1)^2*hypergeom([2-n, 2-n], [2-2*n], 3/4)+((-(1/32)*3^(2+2*n)*n^2-(3/8)*9^n*(n^4+2*n^3-(1/2)*n-1/4))*hypergeom([2+n, 2+n], [2+2*n], 3/4)+9^n*(1/2+n)^2*hypergeom([n+1, n+1], [1+2*n], 3/4)*(3+n)*(n-1/2))*16^(-n)*sin(Pi*n)*n*GAMMA(n)^2)/(Pi^(1/2)*(1/2+n)*(n-1/2)*GAMMA(3/2+n)*n^2*GAMMA(n)*cos(Pi*n)) end proc

seq(simplify(ee(i)), i=2..15);

14, 106, 838, 6802, 56190, 470010, 3968310, 33747490, 288654574, 2480593546, 21400729382, 185239360178, 1607913963614, 13991107041306

evalf(limit(ee(n),n=1));

2.000000000

Download hg_ex.mw

@RaySierra 

Your worksheet's code is missing `+` as the first argument passed to Continue.

If I simply add that missing argument then your worksheet runs ok in my Maple 2019.2.

restart

 

F := proc (q) options operator, arrow; 100*factorial(400)*(1/100)^q*(99/100)^(400-q)/(factorial(400-q)*factorial(q)) end proc

 

F(5.1); evalhf(F(5.1))

15.19780218

Float(undefined)

 

So, F fails under evalhf, since not all the intermediate quantities
get represented as finite values when using hardware double precision.

There are various workarounds which avoid evalhf mode under plot. Here are a few.

 

plot('evalf[Digits](F(q))', q = 0 .. 20)

UseHardwareFloats := false

plot(F, 0 .. 20); UseHardwareFloats := deduced

Since UseHardwareFloats is an environment variable its
assignment inside the proc does not propagate back to
the higher level at which is was called.
F := proc (q) UseHardwareFloats := false; 100*factorial(400)*(1/100)^q*(99/100)^(400-q)/(factorial(400-q)*factorial(q)) end proc

plot(F, 0 .. 20)

UseHardwareFloats

deduced

Download failed_plot_ac.mw

I don't understand what that Matlab colormap is doing, so I'm not sure how to reproduce that exact coloring.

But you could look at this old Post, and adjust the formula for generating the hue component.

Less fast and a bit more clunky are results from complexplot3d, or densityplot (of argument). Eg,

plots:-complexplot3d(ln((exp(1/z))), z=-0.5-0.5*I..0.5+0.5*I,
                     style=surface, orientation=[-90,0,0], grid=[501,501],
                     lightmodel=none, size=[600,600]);
plots:-densityplot(argument(exp(1/(x+I*y))), x=-0.5..0.5, y=-0.5..0.5,
                   axes=boxed, colorbar=false, grid=[301,301],
                   #colorstyle=HUE
                   colorscheme=["zgradient",["Yellow","Blue"],colorspace="HSV"]
                 );

You might compute the derivatives by using the D operator.

When plotted, evalf of such D calls may go through fdiff to perform numeric differentiation. Such approximation of a derivative -- by a finite difference scheme -- of a pde solution numerically approximated via finite difference methods, may need care. I would not expect high accuracy, even if you do try and adjust the step-size(s).

Double-check everything, and adjust as required.

  restart;

  inf:=10:

  pdes:= diff(u(y,t),t)-xi*diff(u(y,t),y)=diff(u(y,t),y$2)/(1+lambda__t)+Gr*theta(y,t)+Gc*C(y,t)-M*u(y,t)-K*u(y,t),
         diff(theta(y,t),t)-xi*diff(theta(y,t),y)=1/Pr*diff(theta(y,t),y$2)+phi*theta(y,t),
         diff(C(y,t),t)-xi*diff(C(y,t),y)=1/Sc*diff(C(y,t),y$2)-delta*C(y,t)+nu*theta(y,t):
  conds:= u(y,0)=0, theta(y,0)=0, C(y,0)=0,
          u(0,t)=0, D[1](theta)(0,t)=-1, D[1](C)(0,t)=-1,
          u(inf,t)=0, theta(inf,t)=0, C(inf,t)=0:
  pars:= { Gr=1, Gc=1, M=1, nu=1, lambda__t=0.5,
           Sc=0.78, delta=0.1, phi=0.5, K=0.5, xi=0.5
         }        

{Gc = 1, Gr = 1, K = .5, M = 1, Sc = .78, delta = .1, nu = 1, phi = .5, xi = .5, lambda__t = .5}

  PrVals:=[0.71, 1.00, 3.00, 7.00]:
  colors:=[red, green, blue, black]:
  for j from 1 to numelems(PrVals) do
      pars1:=`union`( pars, {Pr=PrVals[j]}):
      pdSol:= pdsolve( eval([pdes], pars1),
                       eval([conds], pars1),
                       numeric, timestep=0.01
                     );
      plt[j]:=pdSol:-plot( diff(u(y,t),y), y=0, t=0..2, numpoints=200, color=colors[j]);
  od:

U := (Y,T)->eval(u(y,t),pdSol:-value(t=T)(Y));
U(1,0.5);

Let's check that U produces the same 3D surface as u(y,t)

plots:-display(
  plot3d(U, 0..1, 0..2, style=point, color=orange),
  pdSol:-plot3d(u(y,t), y=0..1, t=0..2, color=blue)
);

proc (Y, T) options operator, arrow; eval(u(y, t), (pdSol:-value(t = T))(Y)) end proc

HFloat(0.0362341204880735)

Plot the derivative of U with respect to its first procedure-parameter (ie, y),
evaluated at y=0.8, for t from 0 to 2.

plot(D[1](U)(0.8,t), t=0..2, color=green)

plots:-display(
  pdSol:-plot3d( u(y,t), y=0..1, t=0..2, numpoints=200, color=colors[1]),
  plottools:-transform((t,dudy)->[0.8,t,dudy])(plot(D[1](U)(0.8,t), t=0..2, color=green)),
  plottools:-transform((t,dudy)->[0.8,t,dudy])(plot(D[2](U)(0.8,t), t=1e-5..2, color=blue)),
  plottools:-transform((y,dudt)->[y,1.8,dudt])(plot(D[2](U)(y,1.8), y=0..1, color=green)),
  plottools:-transform((y,dudt)->[y,1.8,dudt])(plot(D[1](U)(y,1.8), y=0..1, color=blue))
);

## You might play with the space- and time-step, etc.

  PrVals:=[0.71, 1.00, 3.00, 7.00]:
  colors:=[red, green, blue, black]:
  plt:='plt':
  for j from 1 to numelems(PrVals) do
      pars1:=`union`( pars, {Pr=PrVals[j]}):
      pdSol:= pdsolve( eval([pdes], pars1),
                       eval([conds], pars1),
                       numeric, timestep=0.01
                     );
      U := (Y,T)->eval(u(y,t),pdSol:-value(t=T)(Y));
      try
        plt[j]:=plot(D[1](U)(0.00001,t), t=0..2, color=colors[j]);
      catch: end try;
  od:
plots:-display(convert(plt,list));

 

Download badPDE_ac.mw

Try,

   X := Sample(roll,[n,n]):

For Maple 2023 you could use the new colorbar option for 2D contour plots, as well as colorscheme to get that spread of hue values.

For Maple 2022 you could use one of the following approaches:

1) Use DocumentTools:-Tabluate, like in this old response of mine. Adjust weights and dimensions and borders, to taste. You could use the colorscheme option with densityplot, to get the exact colorbar you want, eg. blue to red in hue, etc.

2) Use this old Post of mine, for getting a colored legend along with contourplot (and densityplot too if you want that kind of graded filled effect), or use a specialized form of the legend option with the contourplot command.

3) Do something like the following, using polygons for that filled legend effect. Naturally, you could also merge/display the following with added contour lines.

restart;

with(plots): with(plottools,polygon): with(ColorTools,Color):

 

f := 11-((x-4)^2+(y-3)^2)/25;

11-(1/25)*(x-4)^2-(1/25)*(y-3)^2

(a,b,c,d) := -5,11,-5,11;

-5, 11, -5, 11

Q := plot3d(f, x=a..b, y=c..d):
(zmin,zmax) := (min,max)(op([1,3],Q));

HFloat(5.199999999999999), HFloat(11.0)

display(densityplot(f, x=a..b, y=c..d, 'style'=':-surface',
                        'colorscheme'=["zgradient",["Blue","Red"],
                                       'colorspace'="HSV"]),
  seq(polygon([[0,0],[0,0]],'thickness'=0.3, 'legend'=sprintf("%5.3f",z),
              'color'=Color("HSV",[:-Hue(Color("Blue"))*(zmax-z)/(zmax-zmin),1,1])),
      z=zmax..zmin, 'numelems'=24),
  'legendstyle'=['location'=':-right'], 'size'=[580,500], labels=["",""],
  'axis[2]'=['location'=':-high','tickmarks'=[],'color'="White"],
  'axis[1]'=['location'=':-low','tickmarks'=[],'color'="White"]);

Download M2022_colorbar_fun.mw

I don't know the full details of your example (all ranges, etc), but here's another:

restart;
with(plots): with(plottools,polygon): with(ColorTools,Color):
f := cos(t)*sin(x)-sin(x)+(1/2)*sin(x)*t^2-(1/24)*sin(x)*t^4:
(a,b,c,d) := -Pi,0,-4,1:
Q := plot3d(f, x=a..b, t=c..d):
(zmin,zmax) := (min,max)(op([1,3],Q)):
display(densityplot(f, x=a..b, t=c..d, 'style'=':-surface',
                        'colorscheme'=["zgradient",["Blue","Red"],
                                       'colorspace'="HSV"]),
  seq(polygon([[0,0],[0,0]],'thickness'=0.3, 'legend'=sprintf("%5.3f",z),
              'color'=Color("HSV",[:-Hue(Color("Blue"))*(zmax-z)/(zmax-zmin),1,1])),
      z=zmax..zmin, 'numelems'=24),
  'legendstyle'=['location'=':-right'], 'size'=[580,500], labels=["",""],
  'axis[2]'=['location'=':-high','tickmarks'=[],'color'="White"],
  'axis[1]'=['location'=':-low','tickmarks'=[],'color'="White"]);

I found it by looking at the list of extra options on the ?dsolve,numeric,DAE Help page's Description section.

Near a mention of the differential option, that page had a cross-reference link to the ?dsolve,numeric,DAE_extension Help page, which contains the details of these extension methods for the numeric DAE solvers.

The default extent of the axis view comes from the data points themselves (or a fitted curve).

You can override it by using the usual view option for plotting commands.

For example,

restart; with(Statistics):
N:=200: U:=Sample(Normal(0,1),N):
X,Y:=<seq(1..N)>,<seq(sin(2*Pi*i/N)+U[i]/5,i=1..N)>:
ScatterPlot(X,Y,view=[-20..220,-2..2]);
ScatterPlot(X,Y,view=[-20..220,default]);
ScatterPlot(X,Y,view=[default,-2..2]);

The first sentence in the Options section of the ScatterPlot Help page (Maple 2018) states, with a cross-reference link to the Help page for plot options, "The options argument can contain one or more of the options shown below. All unrecognized options will be passed to the plots[display] command. See plot[options] for details."

What that means is that most of the 2D plotting command common options can be simply passed straight into calls to the ScatterPlot command.

 

restart;

eq:=diff(a(t),t)/a(t)=alpha*t:

T1:=dsolve({eq,a(0)=1}):

de:=diff(lambda(t),t)=rhs(T1):

T2:=dsolve({de,lambda(0)=0}):

SS:=exp(-lambda^2/(2*tau^2))*1/(sqrt(2*Pi)*tau):

SSS:=combine(eval(SS*alpha/rhs(T1),[lambda=rhs(T2),
                                    tau=eval(rhs(T2),t=b)]))
  assuming alpha::positive, b::positive:

expr:=simplify(eval(SSS,[alpha=0.1,b=0.1]));

.3988757910*exp(785.1363895*erf((.2236067978*I)*t)^2-0.5e-1*t^2)

CodeTools:-Usage(
  2*evalf(Int(unapply('evalf[15]'(expr),t),0..infinity,method = _d01amc))
);

memory used=5.23MiB, alloc change=0 bytes, cpu time=37.00ms, real time=36.00ms, gc time=0ns

0.9990021586e-1

Download OverflowError_ac.mw

An alternative to dharr's nice idea, for this example, is to selectively utilize limit.

Naturally, it's not as fast as dharr's switch. But I didn't have to consider the form.

restart;

kernelopts(version);

`Maple 2015.2, X86 64 LINUX, Dec 20 2015, Build ID 1097895`

with(plots):

G0 := y -> exp(y)/(1+exp(y)):
w  := (z, p, q) -> sin(p^+ . z) / (1 + cos(q^+ . z)):

K := 2:
X := Vector(K, symbol=x):
P := `<,>`(2, 1):
Q := `<,>`(0, 3/2):
Modulation := G0(w(X, P, Q)):
pdf        := exp(-1/2 * X^+ . X) / (2*Pi)^(K/2):
Gpdf       := 2*pdf*Modulation:

Gpdf2 := proc(x1,x2) local temp;
  if not [args]::list(numeric) then return 'procname'(args) end if;
  temp := eval(Gpdf,x[1]=x1);
  try eval(temp,x[2]=x2);
  catch: limit(temp,x[2]=x2,':-right');
  end try;
end proc:

opts := grid=[200, 100], contours=[seq(0.01..0.2, 0.02)],
        color=blue, axes=boxed, gridlines=true:
contourplot(Gpdf2(x[1],x[2]), x[1]=-3..3, x[2]=-3..3, opts);

display(
  seq(implicitplot(Gpdf2(x[1],x[2])=level, x[1]=-3..3, x[2]=-3..3,
                   color=blue, gridrefine=3),
      level in eval(contours,[opts])));

 

Download Discontinuous_contours_ac.mw


ps. for implicitplot this use of gridrefine=3 above seems to produce a slightly nicer result that the original gridrefine=4 does there.

Your other set of contour values is also improved:

opts2 := grid=[200, 100], contours=[seq(0.01..0.04, 0.0025)],
        color=blue, axes=boxed, gridlines=true:
contourplot(Gpdf2(x[1],x[2]), x[1]=-3..3, x[2]=-3..3, opts2);

display(
  seq(implicitplot(Gpdf2(x[1],x[2])=level, x[1]=-3..3, x[2]=-3..3,
                   color=blue, gridrefine=3),
      level in eval(contours,[opts2])));

 

First 8 9 10 11 12 13 14 Last Page 10 of 309