## 13653 Reputation

19 years, 286 days

## DirectSearch...

The package DirectSearch is available from the Maple Application Center:
https://www.maplesoft.com/Applications/

```DirectSearch:-GlobalOptima(G(x,y),[x=0..1,y=0..1]); # Minimum
-DirectSearch:-GlobalOptima(-G(x,y),[x=0..1,y=0..1]); # Maximum```

## Using PDEtools:-dchange...

PDEtools:-dchange can be used:

```restart;
e:=Int(1/(3*u^4 + u + 3),u);
res1:=PDEtools:-dchange({u=y*x^(1/3)},e,[y],params={x});
res2:=Int(x/(3*x^2*y^4+x*y+3*x^(2/3)),y); # Result from acer's use of IntegrationTools:-Change
simplify(diff(res1-res2,y)); #0
```

## An animation...

Here is an animation in gamma:

The code:

```restart;
local gamma;
Digits:=15;
#gamma:=0.318;
alpha:=-1;
beta:=1;
delta:=0.1
omega:=1.4;
ics := {v(0) = 0, x(0) = 0};
sys := {diff(v(t), t) = -delta*v(t) - alpha*x(t) - beta*x(t)^3 + gamma*cos(omega*t), diff(x(t), t) = v(t)};
sol := dsolve(sys union ics, numeric, parameters=[gamma],maxfun = 0, abserr = 0.1e-17, relerr = 0.1e-12);
with(plots):
sol(parameters=[0.35]);
odeplot(sol,[t,x(t)],0..800,numpoints=10000,size=[1200,default]);
odeplot(sol,[x(t),v(t)],0..800,numpoints=10000,size=[1200,default]);
odeplot(sol,[x(t),v(t)],2200..3000,numpoints=10000,size=[1200,default]);
Q:=proc(gamma,{range::range:=2200..3000}) if not gamma::realcons then return 'procname'(_passed) end if;
sol(parameters=[gamma]);
plots:-odeplot(sol,[x(t),v(t)],range,numpoints=10000,size=[1200,default]);
end proc;
Q(0.318);
animate(Q,[gamma],gamma=0.318..0.35,frames=100,size=[1200,default]);
```

It makes some sense to me to plot [t,x(t).v(t)] too. This is just as converting the the existing nonautonomous system into an automous one in 3 variables.
The code for Q can just be revised in this way:

```Q:=proc(gamma,{range::range:=2200..3000, scene::list:=[x(t),v(t)]})
if not gamma::realcons then return 'procname'(_passed) end if;
sol(parameters=[gamma]);
plots:-odeplot(sol,scene,range,numpoints=10000,_rest);
end proc;
## and then animate in 3 dimensions:
animate(Q,[gamma,scene=[t,x(t),v(t)]],gamma=0.318..0.35,frames=100);

```

The animation:

Maybe a less confusing version: Narrow the range and make t the third variable instead of the first:

```animate(Q,[gamma,scene=[x(t),v(t),t],range=2900..3000],gamma=0.318..0.35,frames=100);
```

When animating set the FPS (frames per second) to 1.

## MultiSeries:-limit...

Using MultiSeries:-limit gives the same result for both versions:

```restart;
res1:=limit(CylinderU(0,CylinderU(0,x)),x=0);
res2:=CylinderU(0,limit(CylinderU(0,x),x=0));
res1M:=MultiSeries:-limit(CylinderU(0,CylinderU(0,x)),x=0);
res2M:=CylinderU(0,MultiSeries:-limit(CylinderU(0,x),x=0));
```

res1<>res1M, but res2=res2M = res1M= CylinderU(0, 1/2*2^(3/4)*sqrt(Pi)/GAMMA(3/4)).

Further evidence can be obtained from the defining differential equation for CylinderU:

```# y'' - (x^2/4 + a)*y = 0; a = 0 here
ode:=diff(y(x),x,x)-x^2/4*y(x)=0;
dsolve(ode); # Uses Bessel functions
y0:=CylinderU(0,0);
y1:=eval(diff(CylinderU(0,x),x),x=0);
ics:=y(0)=y0,D(y)(0)=y1; # initial conditions
sol:=dsolve({ode,ics});
CU:=unapply(rhs(sol),x); # This is CylinderU(0,x)
limit(CU(CU(x)),x=0);
ev1:=evalf[15](%);
CU(limit(CU(x),x=0));
ev2:=evalf[15](%);
ev_res2:=evalf[15](res2);
```

In fact we can prove the identity of rhs(sol) and CylinderU(0,x):

```B:=convert(CylinderU(0,x),Bessel); # Uses BesselJ only
B2:=convert(rhs(sol),BesselJ);
simplify(B-B2) ; # 0
```

## No need to simplify...

There is really no need to simplify. dsolve converts floats ( i.e. numbers like 0.123 ) into rationals (in this case 123/1000).
That makes the result look complicated.
A simpler version of the solution (sol) is then provided by evalf(sol):

```restart;
ode:=diff(f(x),x)=0.0768*f(x)^(2/3)-0.0102*f(x);# f(1)=59.
sol:=dsolve({ode,f(1)=59});
evalf(sol); # Looks simpler
plot(rhs(sol),x=0..5000);
```

## I believe you are right...

@dharr I think you are quite right:

```restart;
ODE := diff(y(t), t) = -0.000032229*y(t)^2 + 0.065843*y(t) - 15.103;
RHS:=subs(y(t)=y,rhs(ODE));
z1,z2:=solve(RHS=0,y);
plot(RHS,y=z1-100..z2+100);
SOL:=dsolve({ODE,y(0)=z});
plot(eval(rhs(SOL),z=z1),t=0..800); # Constant (roughly)
plot(eval(rhs(SOL),z=z1+0.001),t=0..800);
plots:-animate(plot,[rhs(SOL),t=0..800],z=z1+0.001..z1+10,trace=5);
```

Enlarging the image provided I can see that the curve has to satisfy (t,y) = (0,449), thus this would be the curve:

```plot(eval(rhs(SOL),z=449),t=0..350);
```

To get a direction field you can do

`DEtools:-DEplot(ODE,y(t),t=0..350,[y(0)=449],y=0..2000);`

The image is here:

## convert...

Try this:

`convert~(Res2,units,min);`

Result:
Vector([65*Unit('min'), 70*Unit('min'), 75*Unit('min')])

## Another way...

```restart;
pde:=diff(u(x, t), t, t) + 3 + 2*diff(u(x, t), t) + 4*t + x^2 + x^3/3 + diff(u(x, t), t, x, x) + diff(u(x, t), x, x, x, x) = x*t^2;
ode:=y(x)+diff(y(x),x\$2)+x=sin(x);
##########################
p:=proc(eq::equation) local d,fu,res;
if not has(eq,diff) then return eq end if;
d:=indets(indets(eq,specfunc(diff)),function);
fu:=op(remove(has,d,diff));
res:=selectremove(has,(lhs-rhs)(eq),fu);
res[1]=-res[2]  ## minus put on res[2] (see nm's correction below)
end proc:

p(pde);
p(ode);
```

## Use units=mo...

You could do:

```DateDifference(d1, d2, 'units' = mo);
round(%);
```

## option trace...

Have you tried option trace, as in p:=procedure(....) option trace; local ...; etc end proc.
As an example consider this recent example from MaplePrimes https://mapleprimes.com/questions/238264-Period-Unlimited-Decimal-Fraction

where  I have edited JAMET's procedure and here with option trace:

```periode:=proc(r::rational) option trace; local a,b,c,f,i,k,l,p,q,s:
s:=convert(evalf(abs(r)-trunc(r),50),string):
if  s[1]="." then s:=s[2..length(s)] fi:
a:=numer(simplify(r)):
b:=denom(simplify(r)):
if  igcd(b,10)=1 then ## 1
q:=0:
p:=1:
while (10^(p) mod b)<>1 do
p:=p+1 od:
else
f:=ifactors(b)[2]:
k:=0:l:=0:
for i to nops(f) do
if f[i][1]=2 then k:=f[i][2] fi:
if f[i][1]=5 then l:=f[i][2] fi:
od:
c:=b/(2^k*5^l):
q:=max(k,l):
if c=1 then p:=0 else p:=op(2,periode(1/c)) fi:
fi:
[q,p,[seq(parse(s[k]),k=1+q..q+p)]]
end:```

It works Ok for some examples:

```periode(2/3);
periode(50/7);```

But not for all:

`periode(1007/200035);`

## Is this it?...

```periode:=proc(r::rational) local a,b,c,f,i,k,l,p,q,s:
s:=convert(evalf(abs(r)-trunc(r),50),string):
if  s[1]="." then s:=s[2..length(s)] fi:
a:=numer(simplify(r)):
b:=denom(simplify(r)):
if  igcd(b,10)=1 then ## 1
q:=0:
p:=1:
while (10^(p) mod b)<>1 do
p:=p+1 od:
else
f:=ifactors(b)[2]:
k:=0:l:=0:
for i to nops(f) do
if f[i][1]=2 then k:=f[i][2] fi:
if f[i][1]=5 then l:=f[i][2] fi:
od:
c:=b/(2^k*5^l):
q:=max(k,l):
if c=1 then p:=0 else p:=op(2,periode(1/c)) fi:
fi:
[q,p,[seq(parse(s[k]),k=1+q..q+p)]]
end:
periode(2/3);
periode(50/7);
```

Notice that I changed "couvert"  to "convert".

Note: This fails:

periode(1007/200035);

This appears to happen if q = 0.

To see that q = 0 insert a print statement before the parse statement.

```s:=convert(evalf(abs(1007/200035)-trunc(1007/200035),50),string);
q,p:=0, 1818;
[seq(parse(s[k]),k=1+q..q+p)];
```

## allvalues...

Using allvalues on sol gives two results. odetest produces something not zero though.

```restart;
sol:=y(x) = (exp(RootOf(-sin(x)*tanh(1/2*_Z+1/2*c__1)^2+sin(x)+exp(_Z)))+sin(x))/sin(x);
ode:=diff(y(x),x)-cot(x)*(y(x)^(1/2)-y(x)) = 0;

sol1,sol2:=allvalues(sol);
odetest(sol1,ode,y(x));
odetest(sol2,ode,y(x));
```

## Disabling...

I did it. In order to use it you have to agree to the terms of use.
Wanting to test it I agreed to the terms.  After the testing I took that back. That is possible in the same place.

## Use unapply...

If you want to avoid the weird looks of the definition ot totvol in your second run, use unapply instead of what you got:

`totvol := unapply(Pi*0.4^2*ht^2*ht + 1/3*Pi*0.4^2*ht^2*2/3*ht,ht);`

The output is simple:  totvol := ht -> 0.1955555556*Pi*ht^3

## An illustration...

Illustrating vv's point I tried starting solutions at x=1 for several values of y(1) and plotting the results:

```restart;
de1 := diff(y(x), x) = y(x)/x - 5^(-y(x)/x);
S:=seq(rhs(dsolve({de1,y(1)=k})),k=0..20);
plot([S],x=0..1);
plots:-display(seq(plot(S[k],x=0..1),k=1..21),insequence);# An animation
```

The animation is attached:

 1 2 3 4 5 6 7 Last Page 1 of 160
﻿