Introduction to Matlab

11  Differential Equations.

11.1  Scalar ODE's

The final application is the solution of ordinary differential equations. Suppose we want to find the solution u(t) of the ordinary differential equation

du dt = k u(t) (1-u(t)),    t ³ 0;

(3)

This equation is a simple model for the spread of a disease. Consider a herd of P animals, and let u(t) be the proportion of animals infected by the disease after t days. Thus 1-u(t) is the proportion of uninfected animals. New infections occur when infected and uninfected animals meet. The number of such meetings is proportional to u(t) (1-u(t)). Thus a simple model for the spread of the disease is that the new infections per day, du/dt, is proportional u(t) (1-u(t)). That is we obtain equation (3), where the constant k depends on the density of the animals and the infectiousness of the disease. For definiteness, suppose here P = 10000 animals and k = .2. Initially at time t = 0 there is just one infectious animal, and so we add the initial condition

u(0) = 1/1000.

(4)

We want to find how many individuals will be infected over the next 100 days, that is to plot u(t) from t = 0 to t = 100. To solve this problem we need to define a function to calculate the right hand side of (3), and then use the Matlab program ode45. Use the editor to put the following code into the file `disease.m'.

function  u_dot  =  disease( t , u )      %  Place this code in
%                                         %  file disease.m
%  Evaluates the rhs of the ode
 
k = .2 ;
u_dot  =  k * u * ( 1  - u ) ;

(Although t is not used in calculating u_dot, ode45 requires the function disease to have two input arguments (in the correct order). Use help ode45 for more information.)

Having set up the problem, the differential equation is solved and graphed in just one Matlab one line!

  [ t, u ] = ode45( 'disease', [0,100] , 1/10000  );
  plot( t,u,'x')

The function ode45 also has parameters which can be adjusted to control the accuracy of the computation. Mostly, these are not necessary for just plotting the solution to ordinary problems. However ode45 does not guarantee its accuracy, and it is always a good idea to check our original solution by resolving the problem with higher accuracy requirements. For example, we can use ode45 so that the estimated absolute and relative errors are kept less than 1×10^-9. Also to ensure a very smooth graph, values of the solution are computed every half day. This is done with the commands

  options = odeset('abstol',1.e-9,'reltol', 1.e-9) ;
  [tt,uu] = ode45('disease', [0:.5:100], 1/10000, options);
  plot( t,u,'x', tt,uu,'-' )

We immediately see that the solution computed in the first call to ode45 lies on top of the more accurate solution curve. That is the two solutions coincide to graphical accuracy.

For this simple model it is possible to find an analytic formula for u(t). In fact

u(t) = 1 1+cexp(-kt)      where     c = P-1 .

(5)

We can add this formula to the above graph as a final check.

  t0 =  (0:.5:100)' ;
  c = 10000-1 ; k = .2 ; u0 = 1 ./ ( 1 + c*exp( -k*t0 ) ) ;
  plot( t,u,'x', tt,uu,'-', t0,u0,'--' )

The exact formula gives a curve that is the same as the numerical solution. In fact the difference between the two solutions is at most 2.5×10^-7. So we can be confident in using ode45 to solve the problems below.

Exercise 33

(a.)    Check that (5) is a solution to (3). That is, when we differentiate this formula for u(t) we get the right side of (3). And this formula satisfies (4)

(b.)From the graph of u(t), estimate how many days are needed for half the population to be infected? (Hint: The command grid will make this easier to answer.)

(c.)    Suppose now that after the infection is noted at day t = 0 , an inoculation program is started from day 5, and 200 of the uninfected animals are inoculated per day. After t days, the proportion of animals inoculated per day is 200 (t-5)₊/10000. ((t-5)₊ means max{t-5,0}. This could be coded in Matlab as t5 = max([t-5,0]) for example.) These inoculated animals cannot be infected, an new infections occur when an uninfected, uninoculated animal meets an infected animal. Thus the rate of increase of u is now modelled by the differential equation

du dt     =     k u(t)  ( 1-u(t)-.02 (t-5)+)+;

Now, how many animals are uninfected after 100 days with the inoculation program?

(d.) How many animals would be uninfected if only 100 animals could be inoculated each day from day 5? How many animals would have been saved if 100 animals per day were inoculated from day 0?

This simple model problem also has a formula for the exact solution. However in Matlab it is very easy to get a numerical solution; and this can be done even in complicated problems when there is no simple formulae for the solution. Using Matlab we can simply explore the model by altering the parameters and plotting the solutions!

11.2  Order 2 ODE's

As a second example, consider the ordinary differential equation

d2u dt2 +sin(u(t)) = 0 (6) u(0) = p/4 ,   du dt (0) = 0 (7)

This is a more complicated example because the second derivative of u appears, not just the first derivative. (This is called a second order ODE.) The Matlab ODE solvers cannot use second derivatives directly as they only work with first derivatives. However a sneaky trick can be used. If we keep track of both u(t) and du/dt, then second derivatives of u(t) are just first derivatives of du/dt. More formally, let v(t) denote the value of du/dt. Then ( 6) is equivalent to the two equations

du dt = v(t) ,    and     dv dt = -sin(u(t))

Even more elegantly, we let Y(t) be a vector with the two components

Y(t) = é ê ë Y1(t) Y2(t) ù ú û = é ê ë u(t) v(t) ù ú û ,

and so

dY dt = é ê ë du/dt dv/dt ù ú û = é ê ë v(t) -sin(u(t)) ù ú û = é ê ë Y2(t) -sin(Y1(t)) ù ú û

This is now in the form

dY dt = F(t,Y(t)),

where for any vector Y, F(t,Y) is defined by

F(t,Y) = é ê ë Y2 -sin(Y1) ù ú û

Thus to solve (6) and plot u(t) against t we place the following code in the file F.m

function    Y_dot  =  F( t,Y )
%                                 % Note  Y(2) is the Matlab
%  RHS for the 2nd order ODE      % notation for the second
%                                 % compontent of Y.  It does
%                                 % not mean Y(t) at time t=2
Y_dot  =  [  Y(2) ; -sin( Y(1) ) ] ;

Then we use the command

[t,Y]  = ode45( 'F', [0:.25:20], [pi/4;0] );  % Solve ODE and plot
u  =  Y(:,1)  ;  plot(  t, u,'-' )            %   u(t).

Note that ode45 stores the first and second components of the solution in the first and second columns of the output.

Exercise 34

(a.) From the graph, what is the period of the pendulum?

(b.) Draw a graph of the pendulum's path in the (x,y) plane.

(c.) Often when the pendulum swings through small angles, the equation (6) is approximated by

d2u dt2 +u(t) = 0,

u(0) = p/4 , du dt (0) = 0

How different is the solution to this linear approximation to the exact solution? Are the periods the same?