## 13613 Reputation

19 years, 253 days

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

## A bug in int with assuming...

There seems to be a bug in int here. I changed pi into Pi.
#### NOTE The bug is fixed with Physics update 1717.

```restart;
int(x^2,x=0..x); # Allowed in Maple these days
int(x^2,x=0..N); # Obviously OK
## Now the present case:
eq0:=int((1 - omega)^(k + 1), omega = 0 .. 1) = C*int((1 - sigma*sin(2*Pi*N))^k, N = 0 .. N);
simplify(eq0) assuming k::posint; # Surprising!!
A0:=rhs(eq0/C);
simplify(A0) assuming k::posint; # N

## Now don't allow yourself this double use: N upper limit as well as integration variable:
## The upper limit is now N as above but the variable of integration is x:
##
eq:=int((1 - omega)^(k + 1), omega = 0 .. 1) = C*int((1 - sigma*sin(2*Pi*x))^k, x = 0 .. N);
simplify(eq) assuming k::posint; # The integral still unevaluated
A:=rhs(%/C); # The integral
simplify(A) assuming k::posint,N::integer; # No change
############ 12 concrete k's
L:=[seq](A,k=0..11):
L[1..4]; # Viewing 4 of the elements
L2:=simplify(L) assuming N::integer;
pols:=normal(L2/~N);
degree~(pols,sigma);
```

MaplePrimes24-04-04_int_bug.mw

```A0:=int((1 - sigma*sin(2*Pi*N))^k, N = 0 .. N);
simplify(A0) assuming k::posint; # N The bug
L:=[seq](A0,k=0..11):
L[1..5];
simplify(L) assuming N::integer;  #OK
```

Yet another note: simplify without assumptions returns N too:

`simplify(A0);  # N`

## combine pure and simple...

This works without any explicit reference to trig or exp:

```restart;
combine(exp(sin(a)*cos(b))*exp(cos(a)*sin(b)));
```

Output exp(sin(a + b)).
This works in your simplify case:

`simplify(16^(3/2));`

Procedures can be indexed. Here is a conjured up simple example:

```p:=proc(x) local idx;
if type(procname,'indexed') then
idx:=op(procname);
x,idx
else
x
end if
end proc;
p(4);
p[abc](4);
```

## A way to do it...

Clearly if z is unassigned temp_proc(z, 2); will produce the error you see. It cannot determine if z > 2.
Both of the following work:

```plot('temp_proc(z,2)',z=-2..3);# Notice the unevaluation quotes ' '
plot(rcurry(temp_proc,2),-2..3);# A procedural version: no z anywhere
```

To understand rcurry try rcurry(f,2);
###############

You could change your procedure so that it returns unevaluated if it doesn't receive numerical input:

```restart;
fun := piecewise(x+y > 1, (x+y)^2, x-y);

temp_proc := proc(x, y) if not [x,y]::list(realcons) then return 'procname(_passed)' end if;
local out, ind:

#ind := 9: Not used

if x > y then ind := 1 else ind := 0 end if;

if ind = 1 then out := eval(5*fun, {:-x=x, :-y=y}) else out := eval(-5*fun, {:-x=x, :-y=y}) end if:

return(out);
end proc:

temp_proc(z,2);
temp_proc(1,-2);
plot(temp_proc(z,2),z=-2..3);
```

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