acer

an approach...

Answered: acer 32303

October 20 2023

0 3

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.

>

$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$

>

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])

>

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])

>

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)

>

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)

>

>	#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 # );

>

Download ConflcitInequal_acc.mw

one way...

Answered: acer 32303

October 20 2023

4 0

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);

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

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

Download hg_ex.mw

missing argument...

Answered: acer 32303

October 20 2023

0 0

@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.

Hardware vs software floats...

Answered: acer 32303

October 18 2023

2 1

>

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

>

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.

>

>	Since UseHardwareFloats is an environment variable its assignment inside the proc does not propagate back to the higher level at which is was called.

>

Download failed_plot_ac.mw

colormap...

Answered: acer 32303

October 18 2023

1 4

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"]
                 );

D...

Answered: acer 32303

October 17 2023

0 9

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) );

>	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

Sample(roll,[n,n])...

Answered: acer 32303

October 17 2023

0 4

Try,

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

three ideas for Maple 2022...

Answered: acer 32303

October 16 2023

1 10

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;

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

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

>

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"]);

?dsolve,numeric,DAE_extension...

Answered: acer 32303

October 14 2023

1 3

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.

axis view...

Answered: acer 32303

October 13 2023

0 1

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.

another...

Answered: acer 32303

October 12 2023

2 1

>

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/(2tau^2))1/(sqrt(2Pi)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]));

>	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

Download OverflowError_ac.mw

alternative...

Answered: acer 32303

October 09 2023

2 0

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);

>	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) / (2Pi)^(K/2): Gpdf := 2pdf*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:

indexing...

Answered: acer 32303

October 06 2023

0 0

You can do this simply by indexing into C.

You shouldn't be using the deprecated linalg:-bloackmatrix to construct C. That also produces a deprecated lowercase matrix, by the way. You can use the <...> constructor instead, which here produces a Matrix.

I've edited my response to include the possibility that you'd want to make the V[i] into columns of a Matrix.

>

restart;

>	V[1],V[2],V[3] := Vector([3, 2, -4]), Vector([1, 2/3, 7]), Vector([-9, 0, 1/2]);

V[1], V[2], V[3] := Vector(3, {(1) = 3, (2) = 2, (3) = -4}), Vector(3, {(1) = 1, (2) = 2/3, (3) = 7}), Vector(3, {(1) = -9, (2) = 0, (3) = 1/2})

>	C := <V[1],V[2],V[3]>;

C := Matrix(9, 1, {(1, 1) = 3, (2, 1) = 2, (3, 1) = -4, (4, 1) = 1, (5, 1) = 2/3, (6, 1) = 7, (7, 1) = -9, (8, 1) = 0, (9, 1) = 1/2})

>	C[1..3], C[4..6], C[7..9];

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

>	seq(C[(i-1)3+1..i3], i=1..3);

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

>	F := n->C[(n-1)3+1..n3]:

>	F(1),F(2),F(3);

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

>	M := <V[1]\|V[2]\|V[3]>

M := Matrix(3, 3, {(1, 1) = 3, (1, 2) = 1, (1, 3) = -9, (2, 1) = 2, (2, 2) = 2/3, (2, 3) = 0, (3, 1) = -4, (3, 2) = 7, (3, 3) = 1/2})

>	M[..,1], M[..,2], M[..,3]

Vector(3, {(1) = 3, (2) = 2, (3) = -4}), Vector(3, {(1) = 1, (2) = 2/3, (3) = 7}), Vector(3, {(1) = -9, (2) = 0, (3) = 1/2})

Download Vec_indx.mw

mode=log...

Answered: acer 32303

October 04 2023

0 0

Sure, you could apply log10 to the x and (discarding negatives) the y values. Then force new tickmarks by powering 10, then set the axis locations to `low`, etc. But getting realy nice smooth curves by adaptive sampling becomes tricky. Why reinvent that wheel?

You can compare the loglogplot result with a plot call specifying mode=log for both axes.

I've added gridlines. Notices that the log mode does occur for the vertical axis, though the values are close enough to make the effect less apparent.

loglog_ac.mw

%log...

Answered: acer 32303

October 04 2023

2 3

The expressions log[10](R) and log10(R) will both evaluate to ln(R)/ln(10). So the change occurs even before solve get's its arguments.

You may use the inert version instead. You can utilize the value command later, if you'd like.

For example (using %log[10] since you used log[10], though it's similar with %log10 ),