acer

32303 Reputation

29 Badges

19 years, 308 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

You cannot do it directly as some option to Tabulate.

But you can programmatically get a zoom by using some of the internal commands which Tabulate uses, and making use of an additional option of the InlinePlot constructor.

restart;
with(plots,display): with(plottools,cylinder):

c1 := display(cylinder([1, 1, 1], 1, 3), orientation = [45, 70],
              scaling = constrained, color = red, size = [300, 300]):
c2 := display(cylinder([1, 1, 1], 1, 3), orientation = [45, 70],
              scaling = constrained, color = blue, size = [300, 300]):

with(DocumentTools): with(DocumentTools:-Layout):

InsertContent(Worksheet(
  Table(exterior = none, interior = none,
        Row(Cell(InlinePlot(c1, scale = 1.7)),
            Cell(InlinePlot(c2, scale = 0.4)))))):

However, I'm not sure that the size=[300,300] will be preserved properly upon worksheet Close/Reopen. It's possible the DocumentTools:-Components could also be used, with the Layout:-InlinePlot call further wrapped in a Components:-Plot call that forced the 300x300 dimensions as width and height. Then suppress the component boundaries, etc. (Sorry, I was once more familiar with these subtleties, but it's been a while...)

I would not be surprised if Edgardo puts a fix into a Physics update.

But, here is a hot patch (for Maple 2024.0) that lets ODESteps handle that example.

restart;

kernelopts(opaquemodules=false):
unprotect(Student:-ODEs:-IVPOrder1):
Student:-ODEs:-IVPOrder1:=FromInert(subs(_Inert_FUNCTION(_Inert_ASSIGNEDNAME("indets", "PROC", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_EXPSEQ(_Inert_LOCAL(6), _Inert_NAME("name", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected")))))))=_Inert_FUNCTION(_Inert_ASSIGNEDNAME("indets", "PROC", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_EXPSEQ(_Inert_LOCAL(6), _Inert_FUNCTION(_Inert_NAME("And", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_EXPSEQ(_Inert_NAME("name", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_FUNCTION(_Inert_NAME("Not", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_EXPSEQ(_Inert_NAME("constant", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))))))))),ToInert(eval(Student:-ODEs:-IVPOrder1)))):
protect(Student:-ODEs:-IVPOrder1):
kernelopts(opaquemodules=true):

 

kernelopts(version);

`Maple 2024.0, X86 64 LINUX, Mar 01 2024, Build ID 1794891`

 

ode:=diff(y(x),x)+x*y(x)=1;
ic:=y(0)=0;

diff(y(x), x)+x*y(x) = 1

y(0) = 0

dsolve([ode,ic]);

y(x) = -((1/2)*I)*exp(-(1/2)*x^2)*Pi^(1/2)*2^(1/2)*erf(((1/2)*I)*2^(1/2)*x)

Student:-ODEs:-ODESteps([ode,ic]);

"[[,,"Let's solve"],[,,[(ⅆ)/(ⅆx) y(x)+x y(x)=1,y(0)=0]],["•",,"Highest derivative means the order of the ODE is" 1],[,,(ⅆ)/(ⅆx) y(x)],["•",,"Isolate the derivative"],[,,(ⅆ)/(ⅆx) y(x)=-x y(x)+1],["•",,"Group terms with" y(x) "on the lhs of the ODE and the rest on the rhs of the ODE"],[,,(ⅆ)/(ⅆx) y(x)+x y(x)=1],["•",,"The ODE is linear; multiply by an integrating factor" mu(x)],[,,mu(x) ((ⅆ)/(ⅆx) y(x)+x y(x))=mu(x)],["•",,"Assume the lhs of the ODE is the total derivative" (ⅆ)/(ⅆx) (mu(x) y(x))],[,,mu(x) ((ⅆ)/(ⅆx) y(x)+x y(x))=((ⅆ)/(ⅆx) mu(x)) y(x)+mu(x) ((ⅆ)/(ⅆx) y(x))],["•",,"Isolate" (ⅆ)/(ⅆx) mu(x)],[,,(ⅆ)/(ⅆx) mu(x)=mu(x) x],["•",,"Solve to find the integrating factor"],[,,mu(x)=(e)^((x^2)/2)],["•",,"Integrate both sides with respect to" x],[,,[]=∫mu(x) ⅆx+C1],["•",,"Evaluate the integral on the lhs"],[,,mu(x) y(x)=∫mu(x) ⅆx+C1],["•",,"Solve for" y(x)],[,,y(x)=(∫mu(x) ⅆx+C1)/(mu(x))],["•",,"Substitute" mu(x)=(e)^((x^2)/2)],[,,y(x)=(∫(e)^((x^2)/2) ⅆx+C1)/((e)^((x^2)/2))],["•",,"Evaluate the integrals on the rhs"],[,,y(x)=(-(ⅈ sqrt(Pi) sqrt(2) erf(ⅈ/2 sqrt(2) x))/2+C1)/((e)^((x^2)/2))],["•",,"Simplify"],[,,y(x)=-((ⅈ sqrt(Pi) sqrt(2) erf(ⅈ/2 sqrt(2) x)-2 C1) (e)^(-(x^2)/2))/2],["•",,"Use initial condition" y(0)=0],[,,0=C1],["•",,"Solve for" `c__1`],[,,C1=0],["•",,"Substitute" `c__1`=0 "into general solution and simplify"],[,,y(x)=-ⅈ/2 (e)^(-(x^2)/2) sqrt(Pi) sqrt(2) erf(ⅈ/2 sqrt(2) x)],["•",,"Solution to the IVP"],[,,y(x)=-ⅈ/2 (e)^(-(x^2)/2) sqrt(Pi) sqrt(2) erf(ⅈ/2 sqrt(2) x)]]"

 

 

Download odesteps_fail_nm_May_10_2024_ac.mw

One thing that it shows is that the two results of using alone either simplify or combine coincide, since the set of two identical results has collapsed to a set with only one member.

{simplify,combine}(piecewise(x <= 0, 2*ln(1 - x), 0 < x, ln((x - 1)^2))) assuming (x, real);

{ln((x-1)^2)}

[simplify,combine](piecewise(x <= 0, 2*ln(1 - x), 0 < x, ln((x - 1)^2))) assuming (x, real);

 

[ln((x-1)^2), ln((x-1)^2)]

 

The sentence pertaining to this example in the Maple 2024 Advanced Math What's New page is [emphasis mine],

   "In Maple 2024, simplify and combine are both now better at
    recognizing when two piecewise branch values can be merged into one:"

I think that the main gist of that sentence is that simplify and combine are each separately ("alone") better at dealing with this example.

This is not a bug. Rather, it's a frequently asked question.

The sum command relies on Maple's standard evaluation rules for procedure calls, in which all the arguments are evaluated up front. But before p takes on actual integer values the first argument being passed to sum is actually zero.

In contrast, the add command has special evaluation rules, in which case evaluation of its first argument is delayed until the summing-index name p takes on its actual values.

restart;

e[1] := (t,x,y,z)-> vector(4,[1,0,0,0]):
e[2] := (t,x,y,z)-> vector(4,[0,1,0,-6*z^4*(1-z)^3*(1-x^2)^2*(1-y^2)^3*x]):
e[3] := (t,x,y,z)-> vector(4,[0,0,1,-6*z^4*(1-z)^3*(1-x^2)^3*(1-y^2)^2*y]):
e[4] := (t,x,y,z)-> vector(4,[0,0,0,1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3]):

lie := (j,u,v)-> u[1]*diff(v[j],t)-v[1]*diff(u[j],t)+u[2]*diff(v[j],x)-v[2]*diff(u[j],x)
+u[3]*diff(v[j],y)-v[3]*diff(u[j],y)+u[4]*diff(v[j],z)-v[4]*diff(u[j],z):

G:=(t,x,y,z)->matrix(4,4,[
[1,0,0,0],
[0,1+36*z^8*(1-z)^6*(1-x^2)^4*(1-y^2)^6*x^2,36*z^8*(1-z)^6*(1-x^2)^5*(1-y^2)^5*x*y,-6*z^4*(1-z)^3*(1-x^2)^2*(1-y^2)^3*x*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)],
[0,36*z^8*(1-z)^6*(1-x^2)^5*(1-y^2)^5*x*y,1+36*z^8*(1-z)^6*(1-x^2)^6*(1-y^2)^4*y^2,-6*z^4*(1-z)^3*(1-x^2)^3*(1-y^2)^2*y*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)],
[0,-6*z^4*(1-z)^3*(1-x^2)^2*(1-y^2)^3*x*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3),-6*z^4*(1-z)^3*(1-x^2)^3*(1-y^2)^2*y*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3),(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)^2]
]):
'G(t,x,y,z)'=G(t,x,y,z):

cf:=(k,l,m,t,x,y,z) -> seq(lie(p,e[k](t,x,y,z),e[l](t,x,y,z))*G(t,x,y,z)[p,m],p=1..4):

Warning, (in cf) `p` is implicitly declared local

cf1:=(k,l,m,t,x,y,z) -> sum(lie(p,e[k](t,x,y,z),e[l](t,x,y,z))*G(t,x,y,z)[p,m],p=1..4):

cf1(2,4,2,t,x,y,z);

0

cf2:=(k,l,m,t,x,y,z) -> add(lie(p,e[k](t,x,y,z),e[l](t,x,y,z))*G(t,x,y,z)[p,m],p=1..4):

Warning, (in cf2) `p` is implicitly declared local

cf2(2,4,2,t,x,y,z);

-6*(-6*z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^2*x-6*z^4*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x*(-3*z^2*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3-7*z^3*(z-1)^2*(y^2-1)^3*(x^2-1)^3-2*z^3*(7*z-4)*(z-1)*(y^2-1)^3*(x^2-1)^3)-(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)*(-24*z^3*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x+18*z^4*(1-z)^2*(-x^2+1)^2*(-y^2+1)^3*x))*z^4*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)


This is what gets passed to sum, which does not have special evaluation rules.

lie(p,e[2](t,x,y,z),e[4](t,x,y,z))*G(t,x,y,z)[p,2]

0

cf3:=(k,l,m,t,x,y,z) -> sum('lie'(p,e[k](t,x,y,z),e[l](t,x,y,z))*G(t,x,y,z)[p,m],p=1..4):

cf3(2,4,2,t,x,y,z);

-6*(-6*z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^2*x-6*z^4*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x*(-3*z^2*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3-7*z^3*(z-1)^2*(y^2-1)^3*(x^2-1)^3-2*z^3*(7*z-4)*(z-1)*(y^2-1)^3*(x^2-1)^3)-(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)*(-24*z^3*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x+18*z^4*(1-z)^2*(-x^2+1)^2*(-y^2+1)^3*x))*z^4*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)

 

Download bug_sum_ac.mw

When you use the key option the comparison `<` is done using evalb.

But evalb cannot resolve an inequality involving two items (the norms) which are not both of Maple type numeric.

The abs of the entries of your Vectors are all of Maple type numeric, so the infinity-norm works. But your other norm involves sqrt, and the ensuing radicals are not of type numeric. So comparison via evalb is inadequate.

But those radicals are of type realcons, hence applying evalf (in the key) produces corresponding floats, which are of type numeric.

So it's the very same reason as why you used is as your own check to compare 4th and 5th Vectors.

nb. When the key comparison fails an error is not thrown. But you can see the failure if using the equivalent two-element comparator instead of the key. option -- an error gets thrown,

restart:

V := [seq(LinearAlgebra:-RandomVector(2, generator=1..10), k=1..4)];

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 8, (2) = 7}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 6, (2) = 2})]

Sorting wrt L(+oo) norm

sort(V, key=(t -> norm(t, +infinity)));  # correct

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 6, (2) = 2}), Vector(2, {(1) = 8, (2) = 7})]

map(norm,V);
andmap(type, %, numeric);

[5, 8, 5, 6]

true

Sorting wrt L(2) norm

sort(V, key=(t -> norm(t, 2))); # not correct

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 8, (2) = 7}), Vector(2, {(1) = 6, (2) = 2})]

map(norm,V,2);
andmap(type, %, numeric);

[34^(1/2), 113^(1/2), 41^(1/2), 2*10^(1/2)]

false

evalb(norm(V[4], 2) < norm(V[3], 2)); # neither true nor false!

2*10^(1/2) < 41^(1/2)

is(norm(V[4], 2) < norm(V[3], 2));

true

sort(V, key=(t -> evalf(norm(t, 2)))); # correct

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 6, (2) = 2}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 8, (2) = 7})]

map(evalf@norm,V,2);
andmap(type, %, numeric);

[5.830951895, 10.63014581, 6.403124237, 6.324555320]

true

sort(V, (v1,v2) -> norm(v1,2) < norm(v2,2));

Error, sort: comparison function argument, (v1, v2) -> norm(v1,2) < norm(v2,2), must be a function that always returns true or false

sort(V, (v1,v2) -> evalf(norm(v1,2)) < evalf(norm(v2,2)));

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 6, (2) = 2}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 8, (2) = 7})]

sort(V, (v1,v2) -> is( norm(v1,2) < norm(v2,2) ));

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 6, (2) = 2}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 8, (2) = 7})]

PP := map(norm,V,2);
sort(PP); # done by address
sort(PP, key=(u->u));
sort(PP, key=evalf);
sort(PP, (a,b)->is(a<b));
sort(PP, (a,b)->a<b);

[34^(1/2), 113^(1/2), 41^(1/2), 2*10^(1/2)]

[34^(1/2), 41^(1/2), 113^(1/2), 2*10^(1/2)]

[34^(1/2), 41^(1/2), 113^(1/2), 2*10^(1/2)]

[34^(1/2), 2*10^(1/2), 41^(1/2), 113^(1/2)]

[34^(1/2), 2*10^(1/2), 41^(1/2), 113^(1/2)]

Error, sort: comparison function argument, (a, b) -> a < b, must be a function that always returns true or false

 

Download SortingVectors_ac.mw

You can do that as follows, using assumptions.

Doing it in two separate steps avoids difficulties in the interplay of assuming and the names within the objects.

restart;

with(Student:-MultivariateCalculus):

A := [0, 0, 0];
B := [a, 0, 0];
C := [a, a, 0];
DD := [0, a, 0];
S := [0, 0, h];
d1 := Line(A, C);
d2 := Line(S, DD);

A := [0, 0, 0]

B := [a, 0, 0]

C := [a, a, 0]

DD := [0, a, 0]

S := [0, 0, h]

Student:-MultivariateCalculus:-Line([0, 0, 0], Vector(3, {(1) = a, (2) = a, (3) = 0}), variables = [x, y, z], parameter = t, id = 1)

_m139664740801120

temp := Distance(d1, d2);

abs(a)^2*abs(h)/(2*abs(a*h)^2+abs(a)^4)^(1/2)

simplify(temp) assuming a>0 and h>0;

a*h/(a^2+2*h^2)^(1/2)

Download SMC_simplify_ex.mw

You do not need to use the symbolic option of simplify (which is dodgy in general).

You could also simplify the call directly, if using assume instead of assuming.

restart;

with(Student:-MultivariateCalculus):

assume(a>0);
assume(h>0);

interface(showassumed=0):

A := [0, 0, 0]: B := [a, 0, 0]: C := [a, a, 0]: DD := [0, a, 0]:

S := [0, 0, h]:

d1 := Line(A, C): d2 := Line(S, DD):

simplify(Distance(d1, d2));

a*h/(a^2+2*h^2)^(1/2)

Download SMC_simplify_ex2.mw

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

restart

with(plots)

with(plottools)

display(disk([1, 1], 5, color = red, style = polygon), circle([1, 1], 5, linestyle = dash, color = blue))

NULL

Download how_do_i_change_the_color_of_the_line_around_a_disk_while_using_plottools_ac.mw

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

restart

with(plots)

with(plottools)

setoptions(size = [400, 400])

display(disk([1, 1], 5, color = pink, style = polygon), circle([1, 1], 5, linestyle = dash, color = blue, thickness = 3))

display(disk([1, 1], 5, color = pink, style = polygon), circle([1, 1], 5, color = pink, thickness = 3), circle([1, 1], 5, linestyle = dash, color = blue, thickness = 3))

 

Download disk_edge_2.mw

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

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;

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

 

proc (Gamma, rho) options operator, arrow; 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) end proc

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

Download PP_D_0_ac2b.mw

 

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

.5110796213

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

((81/4)*RootOf(3*_Z^4-_Z^2)^3-(9/2)*RootOf(3*_Z^4-_Z^2))/(-162*RootOf(3*_Z^4-_Z^2)^3+27*RootOf(3*_Z^4-_Z^2))

Float(undefined)

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

((81/4)*RootOf(3*_Z^4-_Z^2)^3-(9/2)*RootOf(3*_Z^4-_Z^2))/(-162*RootOf(3*_Z^4-_Z^2)^3+27*RootOf(3*_Z^4-_Z^2))

Float(undefined)

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

Float(undefined)

# 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<rho<1) there
# are multiple real roots. The generic RootOf doesn't
# *track* the largest real root, or even the largest positive
# real root. So the central-difference formula used for
# numeric differentiation returns nonsense at Gamma=0, which
# is not surprising.

evalf('D[1]'(f)(0.0,0.5));

57734.90192

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

57734.90192

Download signs_derivatves_bivariate_acc2.mw

 

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

.5110796213

# Fortuitously, the largest real root (when Gamma>-1, Gamma<0)
# is part of the "same curve" as the only real root
# for 0<Gamma . So numeric differentiation for D[1](f) will
# behave itself at Gamma=0 and near Gamma=0.
plot(f(Gamma,1/2),Gamma=-1..3,size=[500,200],thickness=3,adaptive=true);

# Fortuitously, the largest real root (when rho>1, rho<3)
# is part of the "same curve" as the only real root
# for -1<rho<1 . So numeric differentiation for D[2](f) will
# behave itself at rho=1 and near rho=1.
plot(f(1/2,rho),rho=-1..2,size=[500,200],thickness=3,adaptive=true);

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

 

Download signs_derivatves_bivariate_acc2b.mw

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.

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

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

but you've instead got,

    sourcefolder:=cat(kernelopts('C:\Users\ahmed\Downloads\adty_v1_0'),"/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:

   sourcefolder:="C:\\\\Users\\ahmed\\Downloads\\adty_v1_0\\wkptest";

   sourcefolder:="C://Users/ahmed/Downloads/adty_v1_0/wkptest";

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

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

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.

@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

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