acer

32333 Reputation

29 Badges

19 years, 325 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

@C_R You can make a first call to simplify with your expression, and trace `simplify/size` to see with what arguments that gets called.

Then you can restart (or forget as needed), and call `simplify/size` directly with the same argument that you ascertain from the first step. You can set printlevel before that.

I don't see how this would be effectively significantly different from what you are suggesting about some new functionality for a "local" printlevel option.

That second step would bypass printing of all the preliminary/final code executed by simplify, `simplify/do`, etc. And the second call could be terminated with a full color to avoid any typesetting details.

First step, from which all we want is to observe the arguments passed to `simplify/size`. (You could also pass the expanded expression, as a related example.)

restart;
trace(`simplify/size`):
simplify((a + b)*(c + d)+ (e + f)*(g + h)):
{--> enter \`simplify/size\`, args = (a+b)*(c+d)+(e+f)*(g+h)
{--> enter \`simplify/size\`, args = (_z1+_z2)*(_z3+_z4)
<-- exit \`simplify/size\` (now in \`simplify/size\`) = (_z1+_z2)*(_z3+_z4)}
{--> enter \`simplify/size\`, args = (_z5+_z6)*(_z7+_z8)
<-- exit \`simplify/size\` (now in \`simplify/size\`) = (_z5+_z6)*(_z7+_z8)}
<-- exit \`simplify/size\` (now in \`simplify/do\`) = (a+b)*(c+d)+(e+f)*(g+h)}

So now we know that we can just pass the expression itself to `simplify/size`, since we've now seen that's what happens when we call simplify itself upon the expression. (For other examples we might see that other arguments are passed at the inner stage of interest.)

Second step,

restart;
printlevel:=100:
`simplify/size`((a + b)*(c + d)+ (e + f)*(g + h)):
printlevel:=1:

Having said all that, I'll mention that -- personally -- I don't feel that printlevel is the most effective tool. I would always rather use the combination showstat/trace/stopat  with each of those being a crucial part of the mix. Others may feel differently, of course.

Here are two ways to get the value of S(t) from the solution returned by dsolve,numeric , where argument t has some numeric value.

You can actually use odeplot itself to get that other plot. Or you can use either way of using S(t) from the dsolve solution.

restart;

with(plots):

b := 0.4: c := 0.1: n := 10^6: p := 0.5:

deS := diff(S(t), t) = -b*S(t)*I0(t);
deI := diff(I0(t), t) = b*S(t)*I0(t) - c*I0(t);
deR := diff(R(t), t) = c*I0(t);

diff(S(t), t) = -.4*S(t)*I0(t)

diff(I0(t), t) = .4*S(t)*I0(t)-.1*I0(t)

diff(R(t), t) = .1*I0(t)

F := dsolve([deS, deI, deR, S(0) = 1 - p, I0(0) = 1/n, R(0) = p],
            [S(t), I0(t), R(t)], numeric,
            method = rkf45, maxfun = 100000):

odeplot(F, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730,
        colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"],
        labels = ["Time (days)", "  Proportion\nof Population "],
        title = "SIR Model with vaccination", size=[500,300]);

odeplot(F, [[t, b*1/c*S(t)*1/n]], t = 0 .. 730, size=[500,300]);

F(100);

[t = 100., S(t) = HFloat(0.46146837378273076), I0(t) = HFloat(0.018483974421123688), R(t) = HFloat(0.5200486517961457)]

eval(S(:-t),F(100));

HFloat(0.46146837378273076)

Reff := t -> b*1/c*eval(S(:-t),F(t))*1/n:
Reff(100);

HFloat(1.845873495130923e-6)

plot(Reff, 0 .. 730, size=[500,300]);


Alternate

F2 := dsolve([deS, deI, deR, S(0) = 1 - p, I0(0) = 1/n, R(0) = p],
            [S(t), I0(t), R(t)], numeric,
            method = rkf45, maxfun = 100000, output=listprocedure):

FP := eval(S(t),F2):

FP(100);

HFloat(0.46146837378273076)

Reff2 := t -> b*1/c*FP(t)*1/n:
Reff2(100);

HFloat(1.845873495130923e-6)

plot(Reff2, 0 .. 730, size=[500,300]);

 

Download Paras31_2_ac.mw

You've forgotten to load LinearAlgebra (or use the full name of Eigenvalues).

You can use Physics:-diff for x(t) as its second argument.

You had a problem in the second set (after restart) with the mapping of eval.

You incorrectly tried to use simplify as some kind of initializer/postprocessor to the Matrix constructor, by making it a final argument. That's not how the Matrix command works. You can wrap the Matrix call in a simplify call (but that might not even be necessary).

In your second set (after restart) you got muddled up trying to assign to A,B,Q.

SI_MODEL_ac.mw

Should the parameter r of periode be r::fraction so that it doesn't run away on integer input? (Or hard-code a trivial result...)

I notice that the inner list returned by periode can vary from that given by the following procedure, eq. for 1/13 .

F := proc(e::rational) local q,r;
  q := NumberTheory:-RepeatingDecimal(e);
  r := RepeatingPart(q);
  [nops(NonRepeatingPart(q)),nops(r),r];
end proc:

In the loop that forms TanSeq, there an unwanted call to tan around sinh(b*i + a).

You need the angles to be the sequence,
   arctan(sinh(a+i*b))

So TanSeq should just be the sequence of sinh(a+i*b) , given that you then apply arctan to each of those.

Those entries sinh(a+i*b) each (already, by definition) represent tan of an angle. You shouldn't be applying tan to them.

ps. You don't need a loop with iterated list concatenation. You can just use seq.

Equal_Incircles_theorem_ac0.mw

Your procedure phi has a conditional test against,
    type(k,even)

But that might not be what you intended. See comparison below,

restart;

type(k,even) assuming k::even;

false

phi := proc(k,x,L)
  if (type(k,even)) then sqrt(2)*sin(Pi*k*x/L)/sqrt(L)
  else sqrt(2)*cos(Pi*k*x/L)/sqrt(L)
  end if;
end proc:

phi(m,x,L) assuming m::even; # not sin!

2^(1/2)*cos(Pi*m*x/L)/L^(1/2)

phi(m,x,L);

2^(1/2)*cos(Pi*m*x/L)/L^(1/2)


Your phi is always returning the cos result when called
with a nonnumeric first argument, ie. a name or an
assumed name. The assumptions about even/oddness
are not taken into account by your phi.

simplify(int(phi(m,x,L)*_h^2/m2*diff(phi(n,x,L),x,x),x=-L/2..L/2))
 assuming m::posint, n::posint, m::even, n::odd;

(2*I)*n^2*((m+n)*(-1)^((1/2)*m-(1/2)*n)+(m-n)*(-1)^((1/2)*m+(1/2)*n))*_h^2*Pi/(m2*L^2*(m^2-n^2))

eval(%,[m=3,n=2]); # compare with next result

-(48/5)*_h^2*Pi/(m2*L^2)

int(phi(3,x,L)*_h^2/m2*diff(phi(2,x,L),x,x),x=-L/2..L/2);

0

is(k,even) assuming k::even;

true


So let's create phi2 which uses is instead of type.

phi2 := proc(k,x,L)
  if (is(k,even)) then sqrt(2)*sin(Pi*k*x/L)/sqrt(L)
  else sqrt(2)*cos(Pi*k*x/L)/sqrt(L)
  end if;
end proc:

phi2(m,x,L) assuming m::even;

2^(1/2)*sin(Pi*m*x/L)/L^(1/2)

simplify(int(phi2(m,x,L)*_h^2/m2*diff(phi2(n,x,L),x,x),x=-L/2..L/2))
 assuming m::posint, n::posint, m::even, n::odd;

0

int(phi2(3,x,L)*_h^2/m2*diff(phi2(2,x,L),x,x),x=-L/2..L/2);

0

 

If you call phi2 without any even/odd kind of assumptions then
there too it returns the cos result.

 

is(k,even);

false

phi2(m,x,L);

2^(1/2)*cos(Pi*m*x/L)/L^(1/2)

simplify(int(phi(m,x,L)*_h^2/m2*diff(phi(n,x,L),x,x),x=-L/2..L/2))
 assuming m::posint, n::posint;
eval(%,[m=3,n=2]);

-2*n^2*((m+n)*sin((1/2)*Pi*(m-n))+sin((1/2)*Pi*(m+n))*(m-n))*_h^2*Pi/(m2*L^2*(m^2-n^2))

-(48/5)*_h^2*Pi/(m2*L^2)

Download test1_3ac.mw

So, having indicated that your phi might not be behaving as you expected, the integration results from trying the four combinations (all with m<>n) with each of m,n being assumed either even/odd might also bear re-examination.

But I'm not sure whether you are expecting results of zero. i.e. when integrating that phi2 variant under similar assumptions.

I suspect that you're asking how to achieve this programmaticaly, entirely in Maple.

You could try thing kind of thing:

with(FileTools:-Text):

WriteFile("f1234.mpl",
          cat(ReadFile("f1.mpl"),"\n",ReadFile("f2.mpl"),"\n",
              ReadFile("f3.mpl"),"\n",ReadFile("f4.mpl")));

You could also add other filler between the contents, if you'd like. You could also use seq and cat, if you have many files to conjoin and they're named in a numeric pattern.

The "\n" is to get a line-break between the end of the last line in a file and the beginning of the first line in the next file.

The student version of Maple has the same computational functionality as the regular version, if I recall correctly. (As far as I know, what differs is price and licensing.)

I'm not sure what you mean by "package" here, unless by that you mean version or licensing.

Here is one way to handle your example, if your goal is to regain the original factored form.

ee := (sqrt(A+B)*x+sqrt(7-K)*y)^2;

((A+B)^(1/2)*x+(7-K)^(1/2)*y)^2

ff := expand(ee)

x^2*A+x^2*B+2*(A+B)^(1/2)*x*(7-K)^(1/2)*y-y^2*K+7*y^2

simplify(Student:-Precalculus:-CompleteSquare(ff));

((A+B)^(1/2)*x+(7-K)^(1/2)*y)^2


Download Ronan_simp_ex.mw

I personally like to see intercepts, and they can also be utilized here.

Also, I think that it's instructive to programmatically construct the equations from the Matrix/Vector data, if you've already shown such.

And this can get rid of implicitplot, which is a computationally inefficient hammer here.

restart

with(plots); with(plottools)


Define the matrix of coefficients and the vector of constants

A := Matrix([[1, 1], [12, 16]]); b := Vector([10, 136])

Matrix(2, 2, {(1, 1) = 1, (1, 2) = 1, (2, 1) = 12, (2, 2) = 16})

Vector[column](%id = 36893628725706251068)


Solve the linear system

sol := LinearAlgebra:-LinearSolve(A, b)

Vector[column](%id = 36893628725706240356)

eqs := Equate(A.`<,>`(x, y), b)

[x+y = 10, 12*x+16*y = 136]

eqs := `~`[isolate](eqs, y)

[y = 10-x, y = 17/2-(3/4)*x]


Compute x-intercepts

X1, X2 := `~`[solve](eval(eqs, y = 0))[]

10, 34/3


Compute y-intercepts

Y1, Y2 := `~`[solve](eval(eqs, x = 0))[]

10, 17/2

display(line([0, Y1], [X1, 0], color = red), line([0, Y2], [X2, 0], color = blue), pointplot([sol], symbol = circle, symbolsize = 12))

display(plot(rhs(eqs[1]), x = min(0, X1) .. max(0, X1), color = red, adaptive = false, numpoints = 2), plot(rhs(eqs[2]), x = min(0, X2) .. max(0, X2), color = blue, adaptive = false, numpoints = 2), pointplot([sol], symbol = circle, symbolsize = 12))

 

Download linear_systems_icepts.mw

I could also have made the above two lines using a single call to the plot command. But I wanted to see them from intercept-to-intercept.

There are so many ways to do all this, including the construction of the equations.

The 3D case can be handled similarly.

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

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