acer

32363 Reputation

29 Badges

19 years, 332 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

The command plottools:-getdata can extract data from these, just like for usual 2D plots of curves.

restart;

plots:-setoptions(size=[600,200]);

pd_numeric:=(D[2,2])(u_numeric)(x,t)+(D[1,1,1,1])(u_numeric)(x,t)-0*h(x,t)=0:
bc_numeric[1]:=u_numeric(0,t)=0:
bc_numeric[2]:=(u_numeric)(1,t)=0:
bc_numeric[3]:=D[1,1](u_numeric)(0,t)=0:
bc_numeric[4]:=D[1](u_numeric)(1,t)=0:
ic_numeric[1]:=u_numeric(x,0)=0.1*x*(x-1)^2:
ic_numeric[2]:=D[2](u_numeric)(x,0)=0:

sol:=pdsolve(pd_numeric, {seq(bc_numeric[i], i=1..4),
             seq(ic_numeric[i],i=1..2)}, u_numeric(x,t), time=t,
             range=0..1, numeric, spacestep=1/2000, timestep=1/2000):

sol:-plot(t=1, labels=[x,u_numeric]);

sol:-value(t=0.5)(0.75);
eval(u_numeric(x, t),%);

[x = .75, t = .5, u_numeric(x, t) = HFloat(-0.0015866001714942412)]

HFloat(-0.0015866001714942412)

UT[0.5] := eval(u_numeric(x, t),
               sol:-value(t=0.5, output=listprocedure)):

UT[0.5](0.75);
UT[0.5](0.25);
UT[0.5](1.0);

HFloat(-0.0015866001714942412)

HFloat(0.003952648473246894)

HFloat(0.0)

# 100 data points
str := time[real]():
sol:-plot(x=0.75, t=0..2, labels=[t,u_numeric], numpoints=100);
time[real]()-str;
op([1,1,2],%%);

5.979

1 .. 100, 1 .. 2

# 4001 data points, it seems
str := time[real]():
sol:-plot(x=0.75, t=0..2, labels=[t,u_numeric]);
time[real]()-str;
op([1,1,2],%%);

6.467

1 .. 4001, 1 .. 2

UXT := unapply('eval(u_numeric(:-x,:-t),sol:-value(:-t=T)(X))',X,T,numeric):

UXT(0.75, 0.5);
UXT(0.25, 0.5);
UXT(1.0, 0.5);

HFloat(-0.0015866001714942412)

HFloat(0.003952648473246894)

HFloat(0.0)

UXT(0.75, 0.2);
UXT(0.25, 0.2);
UXT(1.0, 0.2);

HFloat(-0.004525893648234543)

HFloat(-0.013666119932391071)

HFloat(0.0)

# less efficient, also 100 data points
str := time[real]():
plot(UXT(0.75,t), t=0..2, adaptive=false, numpoints=100,
     color=red, labels=[t,u_numeric]);
time[real]()-str;

6.245

# less efficient, also 4001 data points
str := time[real]():
plot(UXT(0.75,t), t=0..2, adaptive=false, numpoints=4001,
     color=red, labels=[t,u_numeric]);
time[real]()-str;

7.372

# an amusing way, also 100 data points
str := time[real]():
plot(<<seq(0..2,numelems=100)>|Vector[column](op([1,3],sol:-plot3d(t=0..2,x=0..0.75, grid=[2,100]))[2,..])>,
     color=red, labels=[t,u_numeric]);
time[real]()-str;

6.009

# an amusing way, also 4001 data points
str := time[real]():
plot(<<seq(0..2,numelems=4001)>|Vector[column](op([1,3],sol:-plot3d(t=0..2,x=0..0.75, grid=[2,4001]))[2,..])>,
     color=red, labels=[t,u_numeric]);
time[real]()-str;

6.535

Download pds_num_ex1.mw

You've mistakenly indexed things.

For example you have,
   Theta(t)[q]
instead of,
   Theta[q](t)

And similarly you have,
   diff(Theta(t),t)[q]
instead of,
   diff(Theta[q](t),t)

And so on. Note that Theta[1], Theta[2], Theta[3] are all names in Maple. They are indexed names, but can be used for function calls like Theta[1](t), etc.

And the initial conditions are done incorrectly. You can replace with D syntax to express the derivatives evaluated at 0.

I didn't wait for a symbolic solution. But a numeric solution using default method follows.

mrpicky_DE_syntax.mw

Check the edits for correctness.

The command,
   plottools:-getdata(p1)[3]
will return the 2-column Matrix of the data in the p1 plot.

You can then use the ExportMatrix command as one way to export it to a file with Excel-compatible format.

For example,

   ExportMatrix( "foo.xls", plottools:-getdata(p1)[3], target=Excel )

You should use add instead of sum there.

Also, if the Vector is assigned to m then utilize that name, not L.

In the current Maple version you could just call it as add(m).

For example,

m := Vector[column](3, [2, 1, c]):

add(m);

3+c

add(m[k], k=1..3) ;

3+c


Download Vadd.mw

The reason is that m[k] throws an error upon evaluation, where k is just an unassigned name. and m is a Matrix/Vector/Array.

The sum command will get its arguments (such as the reference m[k]) evaluated up front, which is Maple's usual evaluation rules for procedure calls. In contrast, the add command has special evaluation rules which delay the evaluation of first argument m[k] until k actually takes on the numeric values.

By the way, the following magenta error message is actually a URL. In my Maple clicking on it goes to this link, which contains an example and explanation, and the workaround using add. (note: Not all error message URLs go somewhere useful.)

m:=Vector[column](3, [1, 1, 1]);
sum(m[k], k = 1 .. 3) ;

Vector(3, {(1) = 1, (2) = 1, (3) = 1})

Error, bad index into Vector

Download Vsum_err.mw

Do you mean that you want the data in the userinfo message shown in the second attachment shown below?

In Maple 2020 it could be done programmatically as,

kernelopts(version);

`Maple 2020.2, X86 64 LINUX, Nov 11 2020, Build ID 1502365`

eval(`plots/coordRanges`[ellipsoidal],[_a=1,_b=1/2,_c=1/3]);

[[3], [3/4], [1/4], [1 .. 5, 1/2 .. 1, 0 .. 1/2], [0 .. 5, 0 .. 5, 0 .. 5]]


Download coordplot3d_ex1_M2020.mw


And in Maple 2024.0 in can be done as,

restart;

`plots/coordRanges3`[ellipsoidal](1,1/2,1/3);

[[3], [3/4], [1/4], [1 .. 5, 1/2 .. 1, 0 .. 1/2], [0 .. 5, 0 .. 5, 0 .. 5]]

 

infolevel[coordplot3d]:=2:

plots:-coordplot3d(ellipsoidal, size=[300,300]);

plots/coordplot3d: u const values: [3]

plots/coordplot3d: v const values: [3/4]

plots/coordplot3d: w const values: [1/4]

plots/coordplot3d: u range: 1 .. 5

plots/coordplot3d: v range: 1/2 .. 1

plots/coordplot3d: w range: 0 .. 1/2

plots/coordplot3d: view: [0 .. 5 0 .. 5 0 .. 5]

Download coordplot3d_ex1.mw

Issuing,
   showstat(plots:-coordplot3d);
can lead one to,
   showstat(`plots/coordplot3d`);

I'm not sure how much timing performance matters to you. If you only need a final plot then you might choose say 1e-4 or 1e-5 as target accuracy for that, and work backwards if reducing earlier working precision helps with performance. I don't know whether you need high accuracy, or faster timing.

Perhaps the following attachment has some useful bits. The re-usable procedures make it more straightforward to compute for arbitrary m-mu pairs. The plot3d command is invoked with a variable range formula (to respect a condition like mu^2>m, etc).

Your list N has about 50 possible values for mu. The "edge" that the 3d point-plot exhibits is coarse and uneven because the granularity of N doesn't allow all these found points to represent the same degree of closeness to the boundary of the inequality. The points themseleves directly provide only a very rough approximation.

But a style=surfacecontour 3D plot can get you a smoother boundary, with an even closeness for each m value. Naturally you could use another contour value, of your choice of boundary height.

And these procedures could be adjusted. A single dsolve call with two "parameters" might be used. Or an outer procedure might accept just one argument (m value), and then do root-finding (bisection or careful NextZero) for optimal mu value at boundary height. Etc.

Or you might interpolate your grid data programmatically, and then poll arbitarily or (with care for nested precision issues) find the boundary or implicit-plot; I didn't do so, but could it you'd like.

LoopError_acc.mw

Is this what you're trying to do: raise the 2D points in p1 to the vertical height (3D) given by the phi formula?

(This makes the 2D and 3D contour values match exactly, as well as allows the 3D plot to actually contain its data -- stored, and maybe a few other niceties of the 2D command.)

Help_3D_view_ac.mw 

Using your 2D contour plot and its colors,

The following attachment works for me using Maple 18.02, if I first increase the Java Heap size. I made it 2048 instead of default 512.

It takes about 45 seconds for my machine to actually save the animated .gif file (before that it appears empty).

Cuve_1X_Anim_with_Circls_Seq_ac1.mw

I only made and exported one single animated .gif file from the worksheet. I didn't test producing more, in the same worksheet session.

But the Resize did in fact work, and the exported .gif file showed that. 

Here's one way to handle an example like that.

The steps are to select the indices in which you're interested, and then to map table-lookup across those indices, and then take the min of the resulting values.

f := rand(1..100):

T := table([seq([R||i=f(),a||i=f(),b||i=f(),b||i=f()][], i=1..3)]);

table( [( R2 ) = 98, ( a1 ) = 45, ( R3 ) = 38, ( R1 ) = 93, ( a2 ) = 59, ( b3 ) = 96, ( a3 ) = 69, ( b2 ) = 100, ( b1 ) = 6 ] )

 

min(map[2](`?[]`,T,select(hastype,[indices(T)],suffixed(a))));

45

 

 

The steps of that are,

 

select(hastype,[indices(T)],suffixed(a));

[[a1], [a2], [a3]]

map[2](`?[]`,T,%)

[45, 59, 69]

min(%);

45


Download anthr_tab_sel2.mw

It's not clear whether you're asking about how to programmatically determine the version of Maple in which code gets run, or whether you're asking about how to determine the version in which a worksheet was last saved.

I'm guessing it's the former.

The command kernelopts(version) returns a name that contains the version number. That could be converted to a string.

You can chop up that string to take only a particular portion, eg.

str:=convert(kernelopts(version),string);

   str := "Maple 2018.2, X86 64 LINUX, Nov 16 2018, Build ID 1362973"

str[7 .. StringTools:-Search(".",str)-1];

            "2018"

You could use either InertForm:-Display or (undocumented) Typesetting:-Typeset for this.

restart;

 

R := 10:

 

plots:-textplot([600,-Pi*floor(R)*ln(R)-0.01,
                 InertForm:-Display(-Pi*floor(r)*ln(r)),
                 font=[Helvitica,bold,14]],labels=[T," "]);

plots:-textplot([600,-Pi*floor(R)*ln(R)-0.01,
                 Typesetting:-Typeset(-Pi*floor(r)*ln(r)),
                 font=[Helvitica,bold,14]],labels=[T," "]);


Download textplot_ts_ex1.mw

Your given example can be handled by the define command.

restart;

 

define(T,T(a::nonunit(algebraic)+b::nonunit(algebraic))=T(a)+T(b),
         T(a::nonunit(algebraic)*x^n::integer)=a*T(x^n),
         'conditional'(T(a::anything)=a*T(1),
                       _type(a,And(Not(identical(1)),freeof(x)))));

 

T(alpha__1*x^2*y+alpha__2*x^4*t+alpha__3*t*y);

t*alpha__2*T(x^4)+y*alpha__1*T(x^2)+alpha__3*t*y*T(1)


Download define_ex.mw

Is this the kind of effect you're after?

Note that there's no need for any isolation/resubstitution to attain the first target form: collect is enough.

Note a sign difference in the second term of the target form.

The same algsubs result can also be obtained directly from the original expression (expanded), without the intermediate target form.

restart;

expr := 3*G*(`&Delta;&gamma;`*H - `&sigma;y`(`&Delta;&gamma;`) + q)/(-q + `&sigma;y`(`&Delta;&gamma;`))^2;

3*G*(`&Delta;&gamma;`*H-`&sigma;y`(`&Delta;&gamma;`)+q)/(-q+`&sigma;y`(`&Delta;&gamma;`))^2

targ := collect(expr,H,simplify);

3*G*`&Delta;&gamma;`*H/(-q+`&sigma;y`(`&Delta;&gamma;`))^2-3*G/(-q+`&sigma;y`(`&Delta;&gamma;`))

RR := G*`&Delta;&gamma;`=y/3*(q - `&sigma;y`(`&Delta;&gamma;`));

G*`&Delta;&gamma;` = (1/3)*y*(-`&sigma;y`(`&Delta;&gamma;`)+q)

algsubs(RR, targ);

-3*G/(-q+`&sigma;y`(`&Delta;&gamma;`))-H*y/(-q+`&sigma;y`(`&Delta;&gamma;`))

ans := simplify(%);

(H*y+3*G)/(-`&sigma;y`(`&Delta;&gamma;`)+q)

algsubs(RR, expand(expr));

-3*G/(-q+`&sigma;y`(`&Delta;&gamma;`))-H*y/(-q+`&sigma;y`(`&Delta;&gamma;`))

simplify(%);

(H*y+3*G)/(-`&sigma;y`(`&Delta;&gamma;`)+q)

sort(ans, q);

(H*y+3*G)/(q-`&sigma;y`(`&Delta;&gamma;`))


Download ceeeb_ex.mw

[edit] Even though your intermediate target form is not necessary to achieve the algsubs result, it could also be obtained with,

   convert(expr,parfrac,q);
   sort(%,order=plex(H,q));   # optional

For some values of parameter Cl the dsolve command may here return a result containing an inert integral, rather than a fully explicit answer.

In that case you could apply the value command. (This aspect of dsolve is mentioned in a bullet point a section of the ?dsolve,details Help page.)

Eg,

     plot(value(sol(t)), t = 0 .. 48)

dsH_ex.mw

note: There are other ways to deal with it (changing the Int integral to use some specific numeric methods, or doing a numeric dsolve, etc.) Let us know if you'd like to see such.

@jrive First I'll show a form similar to that in nm's Answer, and discuss more on how that is not equivalent to the original even if all parameters are real. (See also nm's pertinent note, "sqrt(a)/b is same as sqrt(a/b^2) assuming b>0", by which he qualifies his answer.)

restart;

sol := (-v + sqrt(-4*a^2*R^2 + v^2))/(2*a*omega*L);

(1/2)*(-v+(-4*R^2*a^2+v^2)^(1/2))/(a*omega*L)

alt := evalindets(expand(sol),`*`,
                  u->`if`(membertype(anything^(1/2),u),sqrt(u^2),u));

-(1/2)*v/(a*omega*L)+(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)


Example for which they're not equal.

eval([sol,alt], [R=1,v=1,a=1,L=1,omega=-1]);

[1/2-(1/2)*(-3)^(1/2), 1/2+(1/2)*(-3)^(1/2)]

combine(simplify( sol-alt )) assuming real, a*omega*L>0;

0

combine(simplify( sol-alt )) assuming real, a*omega*L<0;

(-4*R^2*a^2+v^2)^(1/2)/(a*omega*L)

Download rad_comb_ex.mw

Regarding your comment about bringing a term into a radical, consider this simpler example,

ee := sqrt(x^2)/x;

(x^2)^(1/2)/x

simplify(ee) assuming x>0;

1

simplify(ee) assuming x<0;

-1


If we absorb that x^1 into the radical (by squaring it) and make no other changes then we've mistakenly lost the multiplicative sign of x. If x is positive then it's ok. But if x is negative then we've mistakenly lost the negative sign.

With that in mind, here is a variant on your original query, with a form which retains the signum of a*omega*L outside the radical. Notice the broader equivalence, compared to your original target.

restart;

sol := (-v + sqrt(-4*a^2*R^2 + v^2))/(2*a*omega*L);

(1/2)*(-v+(-4*R^2*a^2+v^2)^(1/2))/(a*omega*L)

alt2 := -v/(2*a*omega*L) + sqrt((-4*R^2*a^2 + v^2)/(a^2*omega^2*L^2))/(2*signum(a*omega*L));

-(1/2)*v/(a*omega*L)+(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)/signum(a*omega*L)

combine(simplify(sol - alt2)) assuming real;

0

alt2 assuming a*omega*L>0;

-(1/2)*v/(a*omega*L)+(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)

alt2 assuming a*omega*L<0;

-(1/2)*v/(a*omega*L)-(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)

Download rad_comb_ex3.mw


Note: If, say, x=-3 then by convention sqrt(x^2) taken to be 3, not -3. That positive answer is known as the principal square root. See also the third paragraph here. Or see notes on the square-root case on principal value or n-th roots.

First 24 25 26 27 28 29 30 Last Page 26 of 336