## 30515 Reputation

18 years, 311 days

## 2D Input operator assignment syntax...

Maple

Caution, certain kinds of earlier input can affect the results from using the 2D Input syntax for operator assignment.

 >

 >
 >

The following now produces a remember-table assignment,
instead of assigning a procedure to name f, even though by default
Typesetting:-Settings(functionassign)
is set to true.

 >

 >

 >

 >
 >

With the previous line of code commented-out the following
line assigns a procedure to name f, as expected.

If you uncomment the previous line, and re-execute the whole
worksheet using !!! from the menubar, then the following will

 >

 >

 >

 >

## integration example...

Maple

There is an interesting preprint here, on a method for integration of some pseudo-elliptic integrals.

It was mentioned in a (sci.math.symbolic) usenet thread, which can be accessed via Google Groups here.

Inside the paper, the author mentions that the following is possible for some cases -- to first convert to RootOf form.

 > restart;
 > kernelopts(version);

 > ig:=((-3+x^2)*(1-6*x^2+x^4)^(-1/4))/(-1+x^2);

 > int(ig,x);

 > infolevel[int] := 2:
 > S:=simplify([allvalues(int(convert(ig,RootOf),x))],size):

Stage1: first-stage indefinite integration

Stage2: second-stage indefinite integration

Norman: enter Risch-Norman integrator

Norman: exit Risch-Norman integrator

int/algrisch/int: Risch/Trager's algorithm for algebraic function

int/algrisch/int: entering integrator at time 9.029

int/algrisch/int: function field has degree 4

int/algrisch/int: computation of an integral basis: start time 9.031

int/algrisch/int: computation of an integral basis: end time 9.038

int/algrisch/int: normalization at infinity: start time 9.039

int/algrisch/int: normalization at infinity: end time 9.048

int/algrisch/int: genus of the function field 3

int/algrisch/int: computation of the algebraic part: start time 9.059

int/algrisch/int: computation of the algebraic part: end time 9.060

int/algrisch/int: computation of the transcendental part: start time 9.063

int/algrisch/transcpar: computing a basis for the residues at time 9.068

int/algrisch/residues: computing a splitting field at time 9.068

int/algrisch/transcpar: basis for the residues computed at time 9.103

int/algrisch/transcpar: dimension is 2

int/algrisch/transcpar: building divisors at time 9.300

int/algrisch/transcpar: testing divisors for principality at time 9.605

int/algrisch/goodprime: searching for a good prime at time 9.606

int/algrisch/goodprime: good prime found at time 9.704

int/algrisch/goodprime: searching for a good prime at time 9.704

int/algrisch/goodprime: good prime found at time 9.762

int/algrisch/areprinc: the divisor is principal: time 10.084

int/algrisch/areprinc: the divisor is principal: time 11.833

int/algrisch/transcpar: divisors proven pincipal at time at time 11.833

int/algrisch/transcpar: generators computed at time 11.834

int/algrisch/transcpar: orders are [1 1]

int/algrisch/transcpar: check that the candidate is an actual antiderivative

int/algrisch/transcpar: the antiderivative is elementary

int/algrisch/transcpar: antiderivative is (1/2)*ln((RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^3*_z^3-RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^2*_z^4+RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^2*_z^2-3*_z^3*RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)+RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)*_z-5*_z^2+1)/((_z+1)*(_z-1)*_z^2))+(1/2)*RootOf(_Z^2+1)*ln((RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^2*RootOf(_Z^2+1)*_z^4+RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^3*_z^3-RootOf(_Z^2+1)*RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^2*_z^2+3*_z^3*RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)-5*RootOf(_Z^2+1)*_z^2-RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)*_z+RootOf(_Z^2+1))/((_z+1)*(_z-1)*_z^2))

int/algrisch/int: computation of the transcendental part: end time 12.040

int/algrisch/int: exiting integrator for algebraic functions at time 12.041

 > S[1];

 > S[2];

 > simplify( diff(S[1],x) - ig ), simplify( diff(S[2],x) - ig );

 >

## 3D plot size option...

Maple 2019

The ideas here are to allow 3D plotting commands such as plot3d to handle a `size` option similarly to how 2D plotting commands do so, and for the plots:-display command to also handle it for 3D plots.

The size denotes the dimensions of the inlined plotting window, and not the relative lengths of the three axes.

I'd be interested in any new problems introduced with this, eg. export, etc.

 > restart;
 > # # Using ToInert/FromInert # # This might go in an initialzation file. # try   __ver:=(u->:-parse(u[7..:-StringTools:-Search(",",u)-1]))(:-sprintf("%s",:-kernelopts(':-version')));   if __ver>=18.0 and __ver<=2019.2 then     __KO:=:-kernelopts(':-opaquemodules'=false);     :-unprotect(:-Plot:-Options:-Verify:-ProcessParameters);     __KK:=ToInert(eval(:-Plot:-Options:-Verify:-ProcessParameters)):     __LL:=[:-op([1,2,1,..],__KK)]:     __NN:=:-nops(:-remove(:-type,[:-op([1,..],__KK)],                           ':-specfunc(:-_Inert_SET)'))           +:-select(:-has,[:-seq([__i,__LL[__i]],                                  __i=1..:-nops(__LL))],                     "size")[1][1];     if :-has(:-op([5,2,2,2,1],__KK),:-_Inert_PARAM(__NN)) then       __KK:=:-subsop([5,2,2,2,1]                      =:-subs([:-_Inert_PARAM(__NN)=:-NULL],                               :-op([5,2,2,2,1],__KK)),__KK);       :-Plot:-Options:-Verify:-ProcessParameters:=       :-FromInert(:-subsop([5,2,2,3,1]                   =:-subs([:-_Inert_STRING("size")=:-NULL],                           :-op([5,2,2,3,1],__KK)),__KK));       :-print("3D size patch done");     else       :-print("3D size patch not appropriate; possibly already done");     end if;   else     :-print(sprintf("3D size patch not appropriate for version %a"),__ver);   end if; catch:   :-print("3D size patch failed"); finally   :-protect(:-Plot:-Options:-Verify:-ProcessParameters);   :-kernelopts(':-opaquemodules'=__KO); end try:

 >
 > P := plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, size=[150,150],             font=[Times,5], labels=["","",""]): P;

 > plots:-display(P, size=[300,300], font=[Times,10]);

 > # # inherited from the contourplot3d (the plot3d is unset). # plots:-display(   plots:-contourplot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1,                        thickness=3, contours=20, size=[800,800]),   plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, color="Gray",          transparency=0.1, style=surface) );

 > # Some options should still act as 2D-plot-specific. # try plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, legend="Q");     print("Not OK"); catch: if StringTools:-FormatMessage(lastexception[2..-1])    ="the legend option is not available for 3-D plots" then print("OK"); else print("Not OK"); error; end if; end try;

 >

If this works fine then it might be a candidate for inclusion in an initialization file, so that it's
automatically available.

## 2D contour plot and legend...

Maple

There have been several posts, over the years, related to visual cues about the values associated with particular 2D contours in a plot.

Some people ask or post about color-bars [1]. Some people ask or post about inlined labelling of the curves [1, 2, 3, 4, 5, 6, 7]. And some post about mouse popup/hover-over functionality [1]., which got added as general new 2D plot annotation functionality in Maple 2017 and is available for the plots:-contourplot command via its contourlabels option.

Another possibility consists of a legend for 2D contour plots, with distinct entries for each contour value. That is not currently available from the plots:-contourplot command as documented. This post is about obtaining such a legend.

Aside from the method used below, a similar effect may be possible (possibly with a little effort) using contour-plotting approaches based on individual plots:-implicitplot calls for each contour level. Eg. using Kitonum's procedure, or an undocumented, alternate internal driver for plots:-contourplot.

Since I like the functionality provided by the contourlabels option I thought that I'd highjack that (and the _HOVERCONTENT plotting substructure that plot-annotations now generate) and get a relatively convenient way to get a color-key via the 2D plotting legend.  This is not supposed to be super-efficient.

Here below are some examples. I hope that it illustrates some useful functionality that could be added to the contourplot command. It can also be used to get a color-key for use with densityplot.

 > restart;
 > contplot:=proc(ee, rng1, rng2)   local clabels, clegend, i, ncrvs, newP, otherdat, others, tcrvs, tempP;   (clegend,others):=selectremove(type,[_rest],identical(:-legend)=anything);   (clabels,others):= selectremove(type,others,identical(:-contourlabels)=anything);   if nops(clegend)>0 then     tempP:=:-plots:-contourplot(ee,rng1,rng2,others[],                                 ':-contourlabels'=rhs(clegend[-1]));     tempP:=subsindets(tempP,'specfunc(:-_HOVERCONTENT)',                       u->`if`(has(u,"null"),NULL,':-LEGEND'(op(u))));     if nops(clabels)>0 then       newP:=plots:-contourplot(ee,rng1,rng2,others[],                               ':-contourlabels'=rhs(clabels[-1]));       tcrvs:=select(type,[op(tempP)],'specfunc(CURVES)');       (ncrvs,otherdat):=selectremove(type,[op(newP)],'specfunc(CURVES)');       return ':-PLOT'(seq(':-CURVES'(op(ncrvs[i]),op(indets(tcrvs[i],'specfunc(:-LEGEND)'))),                           i=1..nops(ncrvs)),                       op(otherdat));     else       return tempP;     end if;   elif nops(clabels)>0 then     return plots:-contourplot(ee,rng1,rng2,others[],                               ':-contourlabels'=rhs(clabels[-1]));   else     return plots:-contourplot(ee,rng1,rng2,others[]);   end if; end proc:
 > contplot(x^2+y^2, x=-2..2, y=-2..2,       coloring=["Yellow","Blue"],       contours = 9,       size=[500,400],       legendstyle = [location = right],       legend=true,       contourlabels=true,       view=[-2.1..2.1,-2.1..2.1] );

 > contplot(x^2+y^2, x=-2..2, y=-2..2,       coloring=["Yellow","Blue"],       contours = 17,       size=[500,400],       legendstyle = [location = right],       legend=['contourvalue',\$("null",7),'contourvalue',\$("null",7),'contourvalue'],       contourlabels=true,       view=[-2.1..2.1,-2.1..2.1] );

 > # Apparently legend items must be unique, to persist on document re-open. contplot(x^2+y^2, x=-2..2, y=-2..2,       coloring=["Yellow","Blue"],       contours = 11,       size=[500,400],       legendstyle = [location = right],       legend=['contourvalue',seq(cat(\$(` `,i)),i=2..5),               'contourvalue',seq(cat(\$(` `,i)),i=6..9),               'contourvalue'],       contourlabels=true,       view=[-2.1..2.1,-2.1..2.1] );

 > contplot(x^2+y^2, x=-2..2, y=-2..2,       coloring=["Green","Red"],       contours = 8,       size=[400,450],       legend=true,       contourlabels=true );

 > contplot(x^2+y^2, x=-2..2, y=-2..2,       coloring=["Yellow","Blue"],       contours = 13,       legend=['contourvalue',\$("null",5),'contourvalue',\$("null",5),'contourvalue'],       contourlabels=true );

 > (low,high,N):=0.1,7.6,23: conts:=[seq(low..high*1.01, (high-low)/(N-1))]: contplot(x^2+y^2, x=-2..2, y=-2..2,       coloring=["Yellow","Blue"],       contours = conts,       legend=['contourvalue',\$("null",floor((N-3)/2)),'contourvalue',\$("null",ceil((N-3)/2)),'contourvalue'],       contourlabels=true );

 > plots:-display(   subsindets(contplot((x^2+y^2)^(1/2), x=-2..2, y=-2..2,                       coloring=["Yellow","Blue"],                       contours = 7,                       filledregions),              specfunc(CURVES),u->NULL),   contplot((x^2+y^2)^(1/2), x=-2..2, y=-2..2,       coloring=["Yellow","Blue"],       contours = 7, #grid=[50,50],       thickness=0,       legendstyle = [location=right],       legend=true),   size=[600,500],   view=[-2.1..2.1,-2.1..2.1] );
 >

 > plots:-display(   contplot(x^2+y^2, x=-2..2, y=-2..2,       coloring=["Yellow","Blue"],       contours = 5,       thickness=0, filledregions),   contplot(x^2+y^2, x=-2..2, y=-2..2,       coloring=["Yellow","Blue"],       contours = 5,       thickness=3,       legendstyle = [location=right],       legend=typeset("<=",contourvalue)),   size=[700,600],   view=[-2.1..2.1,-2.1..2.1] );

 > N:=11: plots:-display(   contplot(sin(x)*y, x=-2*Pi..2*Pi, y=-1..1,       coloring=["Yellow","Blue"],       contours = [seq(-1+(i-1)*(1-(-1))/(N-1),i=1..N)],       thickness=3,       legendstyle = [location=right],       legend=true),    plots:-densityplot(sin(x)*y, x=-2*Pi..2*Pi, y=-1..1,       colorscheme=["zgradient",["Yellow","Blue"],colorspace="RGB"],       grid=[100,100],       style=surface, restricttoranges),    plottools:-line([-2*Pi,-1],[-2*Pi,1],thickness=3,color=white),    plottools:-line([2*Pi,-1],[2*Pi,1],thickness=3,color=white),    plottools:-line([-2*Pi,1],[2*Pi,1],thickness=3,color=white),    plottools:-line([-2*Pi,-1],[2*Pi,-1],thickness=3,color=white),    size=[600,500] );

 > N:=13: plots:-display(   contplot(sin(x)*y, x=-2*Pi..2*Pi, y=-1..1,       coloring=["Yellow","Blue"],       contours = [seq(-1+(i-1)*(1-(-1))/(N-1),i=1..N)],       thickness=6,       legendstyle = [location=right],       legend=['contourvalue',seq(cat(\$(` `,i)),i=2..3),               'contourvalue',seq(cat(\$(` `,i)),i=5..6),               'contourvalue',seq(cat(\$(` `,i)),i=8..9),               'contourvalue',seq(cat(\$(` `,i)),i=11..12),               'contourvalue']),    plots:-densityplot(sin(x)*y, x=-2*Pi..2*Pi, y=-1..1,       colorscheme=["zgradient",["Yellow","Blue"],colorspace="RGB"],       grid=[100,100],       style=surface, restricttoranges),    plottools:-line([-2*Pi,-1],[-2*Pi,1],thickness=6,color=white),    plottools:-line([2*Pi,-1],[2*Pi,1],thickness=6,color=white),    plottools:-line([-2*Pi,1],[2*Pi,1],thickness=6,color=white),    plottools:-line([-2*Pi,-1],[2*Pi,-1],thickness=6,color=white),   size=[600,500] );

 >

## 3D plot of numeric ODE with parameter: p...

Maple

While generating a 3D plot of the solution of an ODE with a parameter, I noticed that better performance could be obtained by calling the plot3d command with a procedure argument, done in a special manner.

I don't recall this being discussed before, so I'll share it. (It it has been discussed, and this is widely known, then I apologize.)

I tweaked the initial conditions of the original ODE system, to obtain a non-trivial solution. I don't think that the particular nature of the solution has a bearing on this note.

 > restart;
 > Digits := 15:
 >

The ODE system has two parameters. One, A, will get a fixed value. The

other, U0, will be used like an independent variable for plotting.

 >
 > eq1:= diff(r[1](t), t, t)+(.3293064114+209.6419478*U0)*(diff(r[1](t), t))       +569.4324330*r[1](t)-0.3123434112e-2*V(t) = -1.547206836*U0^2*q(t): eq2:= 2.03*10^(-8)*(diff(V(t), t))+4.065040650*10^(-11)*V(t)       +0.3123434112e-2*(diff(r[1](t), t)) = 0: eq3 := diff(q(t), t, t)+1047.197552*U0*(q(t)^2-1)*(diff(q(t), t))        +1.096622713*10^6*U0^2*q(t) = -2822.855019*A*(diff(r[1](t), t, t)):
 >
 > ics:=r[1](0)=0,D(r[1])(0)=1e-1,V(0)=0,q(0)=0,D(q)(0)=0:
 >
 > res := dsolve({eq1,eq2,eq3, ics},numeric,output=listprocedure,parameters=[A,U0]):

I will call the procedure returned by dsolve, for evalutions of V(t), as the

dsolve numeric solution-procedure in the discussion below.

 > WV := eval(V(t), res):
 >
 > WV(parameters=[A=1e0]):
 >

The goal is to produce a 3D plot of V(t) as a function of t and the parameter U0.

 >
 > tlo,thi := 0.0, 2.0; U0lo,U0hi := 1e-3, 2e-1;

This is the grid size used for plot3d below. It is nothing special.

 > (m,n) := 51,51;

First, I'll demonstrate that a 3D plot can be produced quickly, by populating a
Matrix for floating-point evaluations of V(t), depending on t in the first
Matrix-dimension and on parameter U0 in the second Matrix-dimansion.

The surfdata command is used. This is similar to how plot3d works.

This  computes reasonably quickly.

But generating the numeric values for U0 and t , based on the i,j positions

in the Matrix, involves the kind of sequence generation formulas that are

error prone for people.

 > str := time[real](): M:=Matrix(m,n,datatype=float[8]): for j from 1 to n do   u0 := U0lo+(j-1)*(U0hi-U0lo)/(n-1);   WV(parameters=[U0=u0]);   for i from 1 to m do     T := tlo+(i-1)*(thi-tlo)/(m-1);     try       M[i,j] := WV(T);     catch:       # mostly maxfun complaint for t above some value.       M[i,j] := Float(undefined);     end try;   end do: end do: plots:-surfdata(M, tlo..thi, U0lo..U0hi,                 labels=["t","U0","V(t)"]); (time[real]()-str)*'seconds';

So let's try it using the plot3d command directly. A 2-parameter procedure
is constructed, to pass to plot3d. It's not too complicated. This procedure
will uses one of its numeric arguments to set the ODE's U0 parameter's

value for the dsolve numeric solution-procedure, and then pass along

the other numeric argument as a t value.

It's much slower than the surfdata call above..

 > VofU0 := proc(T,u0)        WV(parameters=[U0=u0]);        WV(T);      end proc:
 > str := time[real](): plot3d(VofU0, tlo..thi, U0lo..U0hi,        grid=[m,n], labels=["t","U0","V(t)"]); (time[real]()-str)*'seconds';

One reason why the previous attempt is slow is that the plot3d command

is changing values for U0 in its outer loop, and changing values of t in its

inner loop. The consequence is that the value for U0 changes for every

single evaluation of the plotted procedure. This makes the dsolve numeric

solution-procedure work harder, by losing/discarding prior numeric
solution details.

The simple 3D plot below demonstrates that the plot3d command chooses
x-y pairs by letting its second supplied independent variable be the one
that changes in its outer loop. Each time the value for y changes the counter

goes up by one.

 > glob:=0: plot3d(proc(x,y) global glob; glob:=glob+1; end proc,        0..3, 0..7, grid=[3,3],        shading=zhue,  labels=["x","y","glob"]);

So now let's try and be clever and call the plot3d command with the two
independent variables reversed in position (in the call). That will make

the outer loop change t instead of the ODE parameter U0.

We can use the transform command to swap the two indepenent
axes in the plot, if we prefer the axes roles switched. Or we could use the
parametric calling sequence of plot3d for the same effect.

The problem is that this is still much slower!

 > VofU0rev := proc(u0,T)        WV(parameters=[U0=u0]);        WV(T);      end proc:
 > str := time[real](): Prev:=plot3d(VofU0rev, U0lo..U0hi, tlo..thi,              grid=[n,m], labels=["U0","t","V(t)"]): (time[real]()-str)*'seconds'; plots:-display(   plottools:-transform((x,y,z)->[y,x,z])(Prev),   labels=["t","U0","V(t)"],   orientation=[50,70,0]);

There is something else to adjust, to get the quick timing while using

the plot3d command here.

It turns out that setting the parameter's numeric value in the
dsolve numeric solution-procedure causes the loss of previous details
of the numeric solving, even if the parameter's value is the same.

So calling the dsolve numeric solution-procedure to set the parameter

value must be avoided, in the case that the new value is the same as

the old value.

One way to do that is to have the plotted procedure first call the

dsolve numeric solution-procedure to query the current parameter

value, so as to not reset the value if it is not changed. Another way

is to use a local of an appliable module to store the running value

of the parameter, and check against that. I choose the second way.

And plot3d must still be called with the first independent variable-range

as denoting the ODE's parameter (instead of the ODE's independent

variable).

And the resulting plot is fast once more.

 > VofU0module := module()        local ModuleApply, paramloc;        ModuleApply := proc(par,var)          if not (par::numeric and var::numeric) then            return 'procname'(args);          end if;          if paramloc <> par then            paramloc := par;            WV(parameters=[U0=paramloc]);          end if;          WV(var);        end proc: end module:

For fun, this time I'll use the parameter calling sequence to flip the

axes, instead of plots:-transform. That's just because I want t displayed

on the first axis. But for the performance gain, what matters is that it

is U0 which gets values from the first axis plotting-range.

 > str := time[real](): plot3d([y,x,VofU0module(x,y)], x=U0lo..U0hi, y=tlo..thi,        grid=[n,m], labels=["t","U0","V(t)"]); (time[real]()-str)*'seconds';

And, naturally, I could also use the parametric form to get a fast plot

with the axes roles switched.

 > str := time[real](): plot3d([x,y,VofU0module(x,y)], x=U0lo..U0hi, y=tlo..thi,        grid=[n,m], labels=["U0","t","V(t)"]); (time[real]()-str)*'seconds';

 >