Preben Alsholm

13603 Reputation

22 Badges

19 years, 205 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

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

 

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

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

 

Try this:
 

convert~(Res2,units,min);

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

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

 

You could do:

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

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

 

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

 

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

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.

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

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:

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

Note added:
 

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

 

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

 

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);
###############
Note added:

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