Preben Alsholm

13578 Reputation

22 Badges

19 years, 172 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

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

 

You can plot and integrate, but I don't see how you can use dsolve with this object.

restart;
M := Array(1 .. 10, 1 .. 2):

for i to 10 do
    M[i, 1] := i;
    M[i, 2] := 3*i;
end do:

M;
Mt := Interpolation:-LinearInterpolation(M);

Mt(3/2);
plot(Mt(t),t=0..10);
int(Mt,0..10,numeric);

 

The problem is that fun is defined outside the procedure. The x and y there are global variables and has nothing to do with the x and y appearing as the parameters in the procedure.
The remedy is simple:
 

restart;
temp_proc := proc(x, y)
local out1, out2, out3, fun:
fun := x^2+y^2;
if x > 0 then out1 := fun; out2 := 2*fun; out3 := k*fun;
elif x <= 0 then out1 := fun; out2 := -2*fun; out3 := -k*fun;
end if:

return(out1, out2, out3);
end proc:

xt := -1: yt := 2:
out1_fin := temp_proc(xt, yt)[1];
out2_fin := temp_proc(xt, yt)[2];
out3_fin := temp_proc(xt, yt)[3];

In the procedure k is treated as global, so out3_fin depends on whether or not k has been assigned a value outside.

maplemint(temp_proc);


 

I wouldn't use DEplot myself. Here is what I would do:

restart;
with(plots):
sys := {diff(x(t), t) = y(t), diff(y(t), t) = -x(t) - 1/2*y(t)}
ics:={x(0) = 1, y(0) = 0};
sol:=dsolve(sys union ics);
solx,soly:=op(subs(sol,[x(t),y(t)]));
plot([solx,soly,t=0..15],labels=[x,y],color=blue);
animate(plot,[[solx,soly,t=0..T],labels=[x,y],color=red,thickness=2],T=0..20);
###################################
## A nonlinear system. We use numeric solution here.
sys2 := {diff(x(t), t) = y(t), diff(y(t), t) = -x(t) - 1/2*y(t)^2}
ics:={x(0) = 1, y(0) = 0};
res:=dsolve(sys2 union ics,numeric);
odeplot(res,[x(t),y(t)],0..200,numpoints=10000,color="DarkGray");
odeplot(res,[x(t),y(t)],0..200,numpoints=10000,color="SkyBlue",frames=50);

MaplePrimes24-04-01_spiral.mw

 

1 2 3 4 5 6 7 Last Page 1 of 159