## 30152 Reputation

18 years, 232 days

## circle and disk...

You could construct them separately, with the disk having no edge.

You can also compare the effect of having the inner color extend exactly as far as that of the edge (circle).

 >
 >
 >
 >
 >

 >

 >

## Analytic solution?...

If one simply shows the plots generated by your code then it doesn't appear that the numeric solutions for u(y) correspond to the purported "Analytic solution".  So, what is that "Analytic solution" formula supposed to represent? How did you obtain it? Be specific and clear.

Here's an adjustment to get results from your printf attempts. I adjusted the code to keep the dsolve,numeric solutions generated in its loop, so that these could be re-used when constructing the tabular data.

You haven't told us what you mean by "error". Do you mean relative absolute error, or absolute error, or something else? Be specific and clear.

HelpAhsan_ac.mw

## yes...

Yes, there is.

Confusingly, that custom initialization file of preliminary Maple commands is named maple.ini on Windows.

See here for Help with that file, for all of Mac/Windows/Linux.

It's confusing on Windows because -- as you mentioned -- the user's GUI preferences file is named Maple.ini , ie. see here.

These two similarly named files have quite different default locations.

## various parts...

There are several stages to this explanation (which is part of what you asked).

The first is to recall/notice that a generic unspecified RootOf (for a polynomial) doesn't necessarily follow a single continuous curve. Even the indexed RootOf's don't do that. The indexed RootOf's are sorted, and the generic RootOf is taken from the end of those values (so doesn't follow any single one of them throughout). I try to illustrate this in the first attachment below.

The second thing is that when D happens to return unevaluated then evalf of that is a hook into fdiff. But D can actually work symbolically on the f that simply returns the RootOf. And for Gamma=0 the call to D[1] can produce a formula that generates an undefined result. Moreover, there's an important distinction between 'D[1](f)'(1/2,1/2) and D[1](f)(1/2,1/2). The former is an unevaluated call, and applying evalf to it will induce an fdiff attempt. The latter uses the formula that D[1](f)'s symbolic differentiation constructs. (Go back and see that you used the quoted 'D[1]'(f) in your original individual test calls, but in your plot3d attempts you used the unquoted D[1](f) variant. That's why your red and blue surfaces differ in their values along the relevant boundary.)

Also the numeric differentiation across the Gamma=0 and rho=1 boundaries will produce garbage results since the generic RootOf formula switches curve at the boundary on account of the appearance of multiple roots (which affects the principal root via sorting). See the second attachment.
ps. This "switching" affects not just parametrized indexed RootOf's. It also occurs with roots in explicit parametrized radical form. It's not just a Maple thing, but rather it's a Math thing.

The third thing is the issue of how to wade through all these difficulties. You might be able to get around the numeric differentiation issues by using slighty narrow regions. But how much narrower? How to know that features won't get missed? It's not very robust. So my third attachment uses the fact that the greatest root in the multiple root cases (of Gamma<0, Gamma near 0, or rho>1, rho near 1) happens to "match" the surface induced by the single root within Gamma>0,-1<rho<1. So one can compute all four roots (via indexed RootOf's), select the real ones, take the max of those, and get a smooth function across the boundaries Gamma=0, rho=1. So then numeric differentian behaves there. At rho=-1 I just went with a slightly narrower range. This is more time consuming, but the end result seems ok.  See third attachment.

And here they are:

 > restart;
 > plots:-setoptions(size=[500,400]):
 > Eq := -8*(rho + 1)^4*Lambda^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*Lambda^3       - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*Lambda^2       - 4*(rho + 1)*Gamma*(rho^2 - 1)*Lambda + Gamma^2*(rho + 1)*(rho - 1)^2;

 > F := (Gamma,rho) -> RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2);
 >

 > # The generic RootOf in `F` does not track a single # curve across the Gamma=0 boundary, above which there # are multiple real positive roots. # plot(F(Gamma,1/2),Gamma=-1..3,size=[500,200],thickness=3);

 > plot([allvalues(solve(eval(Eq,rho=1/2),Lambda))],Gamma=-1..3);

 > # The generic RootOf in `F` does not track a single # curve across the rho=1 boundary, above which there # are multiple real positive roots. # plot(F(1/2,rho),rho=-1..2,size=[500,200],thickness=3);

 > # For rho>1 there are multiple positive roots. # plot([allvalues(solve(eval(Eq,Gamma=1/2),Lambda))],rho=-2..2);

 > restart;
 > f := (Gamma,rho) -> RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2):
 > evalf(f(1.0,0.5));

 > fDfG := (Gamma,rho) -> fdiff(f, [1], [Gamma,rho]):
 > fDfr := (Gamma,rho) -> fdiff(f, [2], [Gamma,rho]):
 > # This does not call fdiff # `D` can deal with procedure `f`. # See also next. D[1](f)(0,1/2); evalf(%);

 > # See also previous. # NB. this symbolically differentiates # This does not call fdiff. D[1](f): lprint(%); %(0,1/2); evalf(%);

(Gamma, rho) -> (-12*(rho+1)^3*(rho-1)*RootOf(-8*(rho+1)^4*_Z^4+12*(rho+1)^3*
Gamma*(rho-1)*_Z^3-5*(rho+1)^2*(-4/5+Gamma^2*rho^2+2*(-2/5-Gamma^2)*rho+Gamma^2
)*_Z^2-4*(rho+1)*Gamma*(rho^2-1)*_Z+Gamma^2*(rho+1)*(rho-1)^2)^3+5*(rho+1)^2*(2
*Gamma*rho^2-4*rho*Gamma+2*Gamma)*RootOf(-8*(rho+1)^4*_Z^4+12*(rho+1)^3*Gamma*(
rho-1)*_Z^3-5*(rho+1)^2*(-4/5+Gamma^2*rho^2+2*(-2/5-Gamma^2)*rho+Gamma^2)*_Z^2-\
4*(rho+1)*Gamma*(rho^2-1)*_Z+Gamma^2*(rho+1)*(rho-1)^2)^2+4*(rho+1)*(rho^2-1)*
RootOf(-8*(rho+1)^4*_Z^4+12*(rho+1)^3*Gamma*(rho-1)*_Z^3-5*(rho+1)^2*(-4/5+
Gamma^2*rho^2+2*(-2/5-Gamma^2)*rho+Gamma^2)*_Z^2-4*(rho+1)*Gamma*(rho^2-1)*_Z+
Gamma^2*(rho+1)*(rho-1)^2)-2*Gamma*(rho+1)*(rho-1)^2)/(-32*(rho+1)^4*RootOf(-8*
(rho+1)^4*_Z^4+12*(rho+1)^3*Gamma*(rho-1)*_Z^3-5*(rho+1)^2*(-4/5+Gamma^2*rho^2+
2*(-2/5-Gamma^2)*rho+Gamma^2)*_Z^2-4*(rho+1)*Gamma*(rho^2-1)*_Z+Gamma^2*(rho+1)
*(rho-1)^2)^3+36*(rho+1)^3*Gamma*(rho-1)*RootOf(-8*(rho+1)^4*_Z^4+12*(rho+1)^3*
Gamma*(rho-1)*_Z^3-5*(rho+1)^2*(-4/5+Gamma^2*rho^2+2*(-2/5-Gamma^2)*rho+Gamma^2
)*_Z^2-4*(rho+1)*Gamma*(rho^2-1)*_Z+Gamma^2*(rho+1)*(rho-1)^2)^2-10*(rho+1)^2*(
-4/5+Gamma^2*rho^2+2*(-2/5-Gamma^2)*rho+Gamma^2)*RootOf(-8*(rho+1)^4*_Z^4+12*(
rho+1)^3*Gamma*(rho-1)*_Z^3-5*(rho+1)^2*(-4/5+Gamma^2*rho^2+2*(-2/5-Gamma^2)*
rho+Gamma^2)*_Z^2-4*(rho+1)*Gamma*(rho^2-1)*_Z+Gamma^2*(rho+1)*(rho-1)^2)-4*(
rho+1)*Gamma*(rho^2-1))

 > # Same as previous. # This does not call fdiff. evalf(D[1](f)(0,1/2));

 > # This calls fdiff. No symbolic differentiation of # (the body of) f occurs. # Numeric differentiation (w.r.t. Gamma), along Gamma=0, # is problematic because for Gamma<0 (and -1

 > # This calls fdiff # Numeric differentiation along Gamma=0 is problematic. fDfG(0,1/2);

 > restart;
 > f := proc(Gamma,rho) local i;    Digits := 2*Digits;    max(remove(type,simplify(fnormal~([seq(evalf(RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2, index=i)),i=1..4)]),'zero'),nonreal)); end proc:
 > evalf(f(1.0,0.5));

 > # Fortuitously, the largest real root (when Gamma>-1, Gamma<0) # is part of the "same curve" as the only real root # for 0

 > # Fortuitously, the largest real root (when rho>1, rho<3) # is part of the "same curve" as the only real root # for -1

 > # 7sec plot3d(f, 0..10, -1+1e-2..1.0,adaptmesh=false,grid=[23,23],        style=surfacecontour, view=0..1.5);

 > # 47sec plot3d((Gamma,rho)->evalf(D[1](f)(Gamma,rho)),        0..10, -1..1.0, adaptmesh=false, grid=[33,33],        style=surfacecontour, view=-3..0);

 > # 43sec plot3d((Gamma,rho)->evalf(D[2](f)(Gamma,rho)),        0..10, -1.0..1.0, adaptmesh=false, grid=[33,33],        style=surfacecontour, view=-5..1);

 > # 50sec plot3d((Gamma,rho)->evalf(D[1,1](f)(Gamma,rho)),        0..10, -1.0+1e-2..1.0, adaptmesh=false, grid=[33,33],        style=surfacecontour, view=-1..1);

 > # 50sec plot3d((Gamma,rho)->evalf(D[2,2](f)(Gamma,rho)),        0..10, -1.0..1.0, adaptmesh=false, grid=[33,33],        style=surfacecontour, view=-10..300);

 >

## subs...

To address your first query: the subs command does literal replacement. It's not for patterned replacement, e.g. catching any calls to s or any calls to s according to some pattern of its arguments. So, no substitution is done for s(n+1,t) is done if you merely instruct subs to use a formula for s(n,t).

Is this the kind of replacement you wanted for the first query?  shift_ac.mw

Apart from subsindets or evalindets, there is also the applyrule command for replacement of subexpressions by some pattern or rule.

## path syntax...

You gave a link to an old Question. There is a statement that my Answer there:

sourcefolder:=cat(kernelopts('homedir'),"/wkptest");

which does not make sense. What you've got there is not a valid kernelopts call. But it also uses the wrong quotes, ie. neither name nor string quotes.

If you want the source to be located somewhere else (under your Windows home directory, say) then you could try it as one of these:

and so on, presuming that you have actually unpacked it in such a location.

## Explore...

For query 2), if you pass Explore a parameter range with integer end-points, like,
rho = -1 .. 1
then the Slider will only snap to integers -1,0,1.

You could alternatively pass a parameter range with float-point end-points, eg,
rho = -1.0 .. 1.0
to have the Slider take on "more of a continuum" of float values.

## example...

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

## example...

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

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

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

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

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

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

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

 >

## LinearAlgebra...

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

## NumberTheory...

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:```

## tan...

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

## type(k,even)...

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;

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

 > phi(m,x,L);

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;

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

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

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

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;

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

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

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

 > is(k,even);

 > phi2(m,x,L);

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

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.

## FileTools:-Text...

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",