Session 08
Differential Equations and Applications
Maple has a limited ability to find analytical solutions of differential equations. The syntax is dsolve(equation,function,options);
> deq1:=diff(y(t),t)=t+y(t);
> deq1sol:=dsolve(deq1,y(t));
Since you obtain the general solution your answer includes a constant of integration, _C1. Check that this gives a solution.
> u:=rhs(deq1sol);
> simplify(subs(y(t)=u,deq1)); # substitute u for y(t) in the expression deq1
Substitute some values for _C1 (initial values) and plot the solution.
> {seq(subs(_C1=j, u), j=[-2,-1,0,1,2])} # make a set of sequence of expressions u
> plot(%, t=-5..5,-5..5); # plot the set
or :
> plot({seq(subs(_C1=j,u),j=[-2,-1,0,1,2])},t=-5..5,-5..5); # make a set of sequence and plot the set
Now we solve a differential equation with a prescribed initial condition.
> deq2sol:=dsolve({deq1,y(0)=0},y(t));
> plot(rhs(deq2sol),t=-4..4,-8..8);
Sometimes the solutions are given implicitly.
> deq3:=diff(y(t),t)=-t/y(t);
> deq3sol:=dsolve(deq3,y(t));
Sometimes Maple can return explicit results even in these cases by using the explicit option.
> deq3sol:=dsolve(deq3,y(t),explicit);
Note it returns two answers.
Maple has an expanded capability for dealing with differential equations if you load the DEtools package.
> with(DEtools);
> ?DEplot;
> dfieldplot(diff(y(t),t)=cos(t-y(t)),y(t),t=-Pi..Pi,y=-3..3); # plot direction field
dfieldplot(deqns, vars, trange, xrange, yrange, options)
Parameters
deqns - list or set of first order ordinary differential equations
vars - list or set of dependent variables
trange - range of the independent variable
xrange - range of the first dependent variable
yrange - range of the second dependent variable
options- (optional) equations of the form keyword=value
> DEplot(diff(y(t),t)=cos(t-y(t)),y(t),t=-Pi..Pi,{y(0)=1}); # Here we specify an initial condition.
> DEplot(diff(y(t),t)=cos(t-y(t)),y(t),t=-Pi..Pi,{y(0)=1},arrows=NONE);
> DEplot(diff(y(t),t)=cos(t-y(t)),y(t),t=-Pi..Pi,{[0,1],[0,0],[0,-1]}); # initial: [0,1],[0,0],[0,-1]
The above is of the form DEplot(deqns, vars, trange, inits), where inits = initial conditions for solution curves
Here we will solve a system of differential equations, dx/dt=y, dy/dt=-x, for a specified time interval [0,Pi], so that when t=0, x=1 and y=0 (initial conditions). x, y are dependent variables, t is independent variable
> DEplot({diff(y(t),t)=-x(t), diff(x(t),t)=y(t)}, {x(t),y(t)}, t=0..Pi, {[0,1,0]});
The above is of the form DEplot(deqns, vars, trange, inits)
Here we solve a system of nonlinear equations, dx/dt=y, dy/dt=-x+.25*y-y^3. Initial conditions : when t=0, x=0, y=0.1
> DEplot({diff(y(t),t)=-x(t)+.25*y(t)-y(t)^3, diff(x(t),t)=y(t)}, {x(t),y(t)}, t=0..50, {[0,0,.1]}, stepsize=1);
The above is of the form DEplot(deqns, vars, trange, inits, options)
Note the effect of increasing the step size.
> DEplot({diff(y(t),t)=-x(t)+.25*y(t)-y(t)^3, diff(x(t),t)=y(t)}, {x(t),y(t)}, t=0..50, {[0,0,.1]}, stepsize=.1);
This equation exhibits what is called a limit cycle. Here we start a trajectory inside and outside the limit cycle.
> DEplot({diff(y(t),t)=-x(t)+.25*y(t)-y(t)^3, diff(x(t),t)=y(t)}, {x(t),y(t)}, t=0..50, {[0,0,.1],[0,1,0]}, stepsize=.2);
Population dynamics
Let P(t) represent a population at time t. The first and simplest assumption concerning the variation of P(t) is that the rate of change of P is proportional to the current value of P.
> eq1:=diff(P(t),t)=r*P(t);
> sol1:=dsolve({eq1,P(0)=Po},P(t));
> nusol1p:=subs({r=1,Po=10},rhs(sol1));
> nusol1n:=subs({r=-1,Po=10},rhs(sol1));
Note that if r>0 the population grows exponentially and if r<0 it decays exponentially. If the growth rate r =0 then the population stays constant.
To take into account that the growth rate actually depends on the population we now consider a model of the form P'=f(P)P. The first is the so-called Logistic model which behaves similarly to the above model for P(t) very small.
> eq2:=diff(P(t),t)=r*(1-P(t)/K)*P(t);
> dsolve({eq2,P(0)=Po},P(t));
We will consider two cases by way of example.
> DEplot(subs({r=1,K=2},eq2),P(t),t=0..3,{[0,.1],[0,1],[0,3]});
Note that for initial condition P(0)=0 we have P(t)=0 (try substitute for [0,0] as initial condition)and if P(0)=K, then P(t)=K for all t. (Try substitute in this case K=2 for the initial condition, i.e., [0,2]) On the other hand for any positive P(0) the solution goes to K (called the saturation level). For this reason P(t)=K is called an asymptotically stable equilibrium and N(t)=0 is called an unstable equilibrium.
> DEplot(subs({r=1,K=2},eq2),P(t),t=0..3,{[0,0],[0,2],[0,.1],[0,1],[0,3]}); # K=2
> DEplot(subs({r=1,K=1},eq2),P(t),t=0..3,{[0,0],[0,1],[0,.1],[0,2],[0,3]}); # K=1
A great many physical problems can be formulated in terms of second order linear differential equations in the form
 ,
u(0)=u0, u'(0)=u1 (initial conditions)
Two such examples are the motion of a mass on a spring and the flow of electric current in a simple series circuit. In the case of a mass m attached to a spring with spring constant k and if we assume a resistance to motion due to the surrounding medium that is proportional to the speed of the mass, then
To obtain a unique solution the equations need to be supplemented by an initial displacement u(0)=u0 and an initial velocity u'(0)=u1.
First consider the case of simple harmonic motion, i.e., no damping
Let omega =sqrt(k/m). Omega determines the frequency of oscillation.
> u:='u';
> eq1:=diff(u(t),t$2)+omega^2*u(t)=0; 
> sol1:=dsolve({eq1,u(0)=u0,D(u)(0)=u1},u(t)); 
> sol:=subs({omega=4,u0=1,u1=1},rhs(sol1)); 
> plot(sol,t=0..Pi);
Thus we see that the mass simply oscillates periodically forever.
Now we assume there is damping due to friction b=r/m which is not zero. It is convenient to write (without loss of generality) the damping as 2*frequency*b.
Initial conditions: u(0)=u0, u'(0)=u1
> dampeq:=diff(u(t),t$2)+2*b*omega*diff(u(t),t)+omega^2*u(t)=0; 
> dsolve({dampeq,u(0)=u0,D(u)(0)=u1},u(t));
There are three cases to consider
case (a). b^2-1=0
case (b). b^2-1>0 and
case (c). b^2-1<0
Case (a) b^2-1=0, b=1
> damsol:=rhs(%);
> damsola:=limit(damsol,b=1); 
Here is a particular case
> psola:=subs({u0=.1,u1=-1,omega=sqrt(4)},damsola); 
> plot(psola,t=0..5);
This case is known as critical damping. The solution decays to rest.
Case (b) b^2-1 <0
Consider case (b) when b^2-1 <0. Choose as example: (b^2=1/50, u0=1, u1=-1) as a particular case.
> psolb:=subs({u0=1,u1=-1,omega=2,b=sqrt(1/50)},damsol);
> plot(psolb,t=0..3*Pi);
Next we plot the phase portrait, i.e., (u(t),u'(t)).
> plot([psolb,diff(psolb,t),t=0..10]);
In this case we obtain a damped oscillation, i.e., the solution oscillates but with decreasing amplitude due to friction.
Case (c) b^2-1>0
The final case when b^2-1>0 is called over-damping. Recall that we are concerned with the equation.
> dampeqd:=subs({omega=2,b=101/100},dampeq); 
> psold:=dsolve({dampeqd,u(0)=2,D(u)(0)=1},u(t));
> plot(rhs(psold),t=0.. 3);
Note that in the case of over damping the solution simply decays to zero with no oscillations.
Exercises
Problem 1
Use dsolve to find the general solution to the differential equation. Also substitute the solution into the equation and show that it is a solution (i.e., check your answer). If there is an initial condition also find the unique solution with the given initial condition.
a. diff(y(x),x)-6*y(x)=10*sin(2*x); with y(0)=0;
b. diff(y(x),x)-(2/x)*y(x)=((x-2)/x)*exp(x); with y(2)=0
c. diff(r(theta),theta)-(2*r(theta)*cot(theta)+sin(2*theta))=0;
Problem 2
Use DEtools sketch the direction fields (dfieldplot) of the first order equation then use DEplot to plot a few trajectories with the given initial conditions and for t from 0 to 2 (note you may need to restrict the y range);
a. diff(y(t),t)=1/2*y(t)*(y(t)-1); [0,.5], [0,1.25] for t from 0 to 2
b. diff(y(t),t)=sin(4*t*y(t)); [0,.1],[0,.5],[0,1.25],[0,1.6] for t from 0 to 2
c. diff(y(t),t)=-sin(x)-y/10, use [t,x,y] as variable and initial condition [0,1,1] for t from -10 to 10 in DEplot
<< Back - Next >>
|
