Preben Alsholm

MaplePrimes Activity


These are answers submitted by Preben Alsholm

You could use odeplot from the plots package.

Here is an example which I have put together hastily (I admit):
 

restart;
sys := {diff(x(t), t) = x(t) - y(t), diff(y(t), t) = 2*x(t) - y(t) - x(t)^2};
DEtools:-DEplot(sys,[x(t),y(t)],t=-2..5,x=-1..3,y=-2..3,[[x(0)=2,y(0)=3],[x(0)=1/2,y(0)=1/2]],linecolor=blue);
ics:={x(0)=2,y(0)=3};
res:=dsolve(sys union ics,numeric);
plots:-odeplot(res,[t,sqrt(x(t)^2+y(t)^2)],0..2);
res(initial);
res(initial=[0,1/2,1]);
plots:-odeplot(res,[t,sqrt(x(t)^2+y(t)^2)],0..8);
?dsolve,numeric,interactive
PL:=proc(IC::listlist,{scene::list:=[x(t),y(t)],range::range:=0..2}) local ic,p;
   for ic in IC do
      res('initial'=ic);
      p[ic]:=plots:-odeplot(res,scene,range,_rest)
   end do;
   plots:-display(seq(p[ic],ic in IC))
end proc;
###Examples:
PL([seq([0,i,1/2],i=-1..1,0.1)]);
PL([seq([0,i,1/2],i=-1..1,0.1)],range=-2.3..2);
PL([seq([0,i,1/2],i=-1..1,0.1)],range=-3..4,view=[-1..3,-2..3]);
PL([seq([0,i,1/2],i=-1..1,0.1)],range=-3..4,scene=[t,sqrt(x(t)^2+y(t)^2)],view=[-3..4,0..10]);

The last plot:

It seems that algsubs won't touch attempts to replace functions like int, Int, sin, or f in f(x):
 

restart;
expr:=int(f(x),x);
algsubs(f(x)=g(x),expr); # Works
algsubs(f=g,expr); # Doesn't do anything
algsubs(sin=cos,sin(x)); # Doesn't do anything
algsubs(x=y,expr); # Works

 

That problem with y(0) = 0 has infinitely many solutions (a well known fact):
 

restart;
ode := diff(y(x),x) = 2*sqrt(y(x));
p:=seq(piecewise(x<n,0,(x-n)^2),n=0..4);
plot([p],x=-1..5);
seq(odetest(y(x)=q,ode),q=p);

With 17 equations we may assume that your solution is numerical.

I give here a very simple example: Just one ode and it can be solved symbolically as well as numerically:
 

restart;
ode1:=diff(x(t),t)=-x(t)+sin(t);
ic1:=x(0)=1;
res1:=dsolve({ode1,ic1},numeric); # The numerical solution for x
plots:-odeplot(res1,[t,x(t)],0..20); 
### Now introduce the integral of x which is called y: 
INT:=y(t)=int(x(s),s=0..t);
### Differentiate that to get an extra ode:
ode2:=diff(INT,t);
ic2:=eval(INT,t=0); # The initial value for y
res2:=dsolve({ode1,ode2,ic1,ic2},numeric); # Numerical solutions for x and y
plots:-odeplot(res2,[[t,x(t)],[t,y(t)]],0..20); 
### Now for comparison we solve exactly (symbolically):
sol1:=dsolve({ode1,ic1});
map(int,sol1,t=0..T);
sol2:=y(t)=subs(T=t,rhs(%));
plot([rhs(sol1),rhs(sol2)],t=0..20); 
### The difference between the numerical and the exact y(t):
plots:-odeplot(res2,[t,y(t)-rhs(sol2)],0..20);

PS. I see that Carl Love beat me to it, but I shall leave my version. It doesn't hurt with several examples especially when you are new to Maple.

Could you not use events?
Here is a suggestion:
 

sysE := {
         diff(x(t), t) = v(t),
         diff(v(t), t) = t - (b(t)*f+(1-b(t))*(-f)),
         x(0) = 0,
         v(0) = 0,
         b(0) = 0
       };
solE := dsolve(sysE, numeric, parameters=[f],discrete_variables=[b(t)::boolean], events=[[v(t)=0,b(t)=1-b(t)]]);
solE(parameters=[1]):
solE(1);
odeplot(solE,[t,v(t)],0..100);

 

In the revised version of U it pays to reverse the order of integration:

 

restart;
#HAM := [1, 2, 3, 4, 5, 6]; 
U := q^6*sin(p); 
alpha := 1/3;
## Notice that the order of integration has been reversed:
for i from 1 to 6 do
   rho := i;
   f[i] := int(int(U/((x^rho-p^rho)^alpha*(y^rho-q^rho)^alpha), p = 0 .. x), q = 0 .. y) assuming x>0,y>0
end do;

 


 

restart;

sys:=P*diff(f(eta), eta, eta, eta)/Q + lambda*theta(eta)*(S0 + S1*Qc*theta(eta)) + (n + 1)/2*(f(eta) + g(eta))*diff(f(eta), eta, eta) - n*diff(f(eta), eta)*(diff(f(eta), eta) + diff(g(eta), eta)) = 0, P*diff(g(eta), eta, eta, eta)/Q + (n + 1)/2*(f(eta) + g(eta))*diff(g(eta), eta, eta) - n*diff(g(eta), eta)*(diff(f(eta), eta) + diff(g(eta), eta)) = 0, (S2 - 2*R)*diff(theta(eta), eta, eta)/(S4*Pr) + (3*R)/(2*S4*Pr*(thetaw - 1))*diff((1 + (thetaw - 1)*theta(eta))^2, eta, eta) + (n + 1)/2*(f(eta) + g(eta))*diff(theta(eta), eta) = 0, f(0) = 0, g(0) = 0, theta(0) = 1, D(f)(0) = 1, D(g)(0) = c, D(f)(3) = 0, D(g)(3) = 0, theta(3) = 0;
 

sys[3]; # Term 2 has a removable singularity at thetaw=1

sys3:=map(normal,lhs(sys[3]))=0; # No singularity at thetaw=1

isolate(sys3,diff(theta(eta),eta,eta));

When R > 0, theta(eta) > 0, and thetaw > 1 we see that the denominator is positive (in fact grater than R+S2).

SYS:=sys[1..2],sys3,sys[4..-1];

P := 0.7952865489: Q := 1.239340734: S4 := 1.007143462:S0 := 0.8330195568:S1 := 0.8330195568:
S2 := 1.087762421: Pr := 0.0157*2415/0.252: n:=3: lambda:=20: Qc:=2: c:=0.5:

for s in SYS do s end do; # Viewing the system

Digits:=15:

Boundary value problems of this kind are notoriously difficult, so we begin with simply providing simple values for thetaw and R:

res:=dsolve(eval({SYS},{thetaw=1.5,R=0.2}),numeric,initmesh=256,maxmesh=1024);

res(0);

dth0:=subs(res(0),diff(theta(eta),eta));

plots:-odeplot(res,[[eta,f(eta)],[eta,g(eta)],[eta,theta(eta)]],0..3);

This procedure will be used for animation:

resplot:=proc(thetaw0,R0) local res;
   res:=dsolve(eval({SYS},{thetaw=thetaw0,R=R0}),numeric,_rest);
   plots:-odeplot(res,[[eta,f(eta)],[eta,g(eta)],[eta,theta(eta)]],0..3,legend=[f,g,theta])
end proc;

resplot(1.5,0.2,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta)]);

plots:-animate(resplot,[1.5,R,initmesh=512,maxmesh=2048,
       approxsoln=[f(eta)=0.4*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta)]],R=0.1..2,trace=24);

plots:-animate(resplot,[thetaw,0.2,initmesh=512,maxmesh=2048,
              approxsoln=[f(eta)=0.4*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta)]],thetaw=1..10,trace=24);

We notice that the general shape of the solutions doesn't change much with R and thetaw.

 

The following procedure takes values of thetaw and D(theta)(0) as input and outputs the corresponding value for R.

Notice the extra boundary condition D(theta)(0)=Dth0:

Rproc:=proc(thetaw0,Dth0) local res;
   if not [thetaw0,Dth0]::list(realcons) then return 'procname(_passed)' end if;
   res:=dsolve(eval({SYS,D(theta)(0)=Dth0},thetaw=thetaw0),numeric,_rest);
   subs(res(0),R)
end proc;

Here we should expect to see the result 0.2 for R even though the R value given in approxsoln is R = 1, and indeed we do:

Rproc(1.5,dth0,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]);

Rproc(1.5,dth0,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=2]);

Rproc(1.5,-1,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]);

plot(
     Rproc(thetaw,-1,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]),
     thetaw=1..2,labels=[thetaw,R]);

 


Here as we saw we have D(theta)(0) = -1 (the second argument to Rproc):

This code:
 

plot3d(Rproc(thetaw,Dth0,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]),
     thetaw=1..2,Dth0=-3..-1,labels=[thetaw,Dth0,R]);

produces (after a while):