AMTH247 Lecture 19

Differential Equations I

Reading: Heath §9.1, 9.2, 9.3.1

Contents

• Theory
  • First Order Equations
  • First Order Systems
  • Higher Order Equations
  
  
• Euler's Method
  • Systems of Equations
  
  
• About this document ...

Theory

First Order Equations

A first order differential equation is an equation of the form

$$\frac{dy}{dt} = f(t, y) .$$ (1)

Once an initial value

$$y(t_0) = y_0$$

is specified, we have an initial value problem. An initial value problem has a unique solution $y(t)$ giving $y$ as a function of $t$. Another way of saying the same thing is that for each point $(t_0, y_0)$ of the $t$-$y$ plane there is a unique solution of equation (1) passing through that point.

Geometrically, the differential equation can be interpreted as saying that the solution curve through a point $(t,y)$ has slope $f(t,y)$ at that point.

First Order Systems

A system of first order differential equations has the form

$$\frac{dy_1}{dt} = f_1 (t, y_1, \dots , y_n)$$

$$\frac{dy_2}{dt} = f_2 (t, y_1, \dots , y_n)$$

$$\vdots$$

$$\frac{dy_n}{dt} = f_n (t, y_1, \dots , y_n)$$

that is, we have a system of $n$ first order differential equations for $n$ functions $y_1(t), \dots, y_n(t)$. An initial value problem specifies the initial values of all $n$ functions at the same point $t_0$:

$$y_1 (t_0) = y_{1,0}$$

$$y_2 (t_0) = y_{2,0}$$

$$\vdots$$

$$y_n (t_0) = y_{n,0}$$

As in the case of a single equation, an initial value problem for a first order system has a unique solution. It is often convenient to write the equations in vector form; if $\mathbf{y} = (y_1, \dots, y_n)$ then the system of differential equations can be written

$$\frac{d\mathbf{y}}{dt} = \mathbf{f}(t, \mathbf{y})$$

where $\mathbf{f}$ is the vector of functions $f_1, \dots, f_n$. The initial values can also be written in vector form

$$\mathbf{y}(t_0) = \mathbf{y}_0 .$$

Higher Order Equations

The differential equations we have studied so far contain only first derivatives. Higher order differential equations involve the second and possibly higher derivatives. An $\mathbf{n}$th-order differential equation has the form

$$\frac{d^ny}{dt^n} = f\left(t, y, \frac{dy}{dt}, \frac{d^2y}{dt^2}, \dots, \frac{d^{n-1}y}{dt^{n-1}} \right) .$$

An initial value problem specifies the initial values of all the derivatives up to order $n-1$:

$$y(t_0) = y_0, \qquad \frac{dy}{dt}(t_0) = {y_0}', \quad \dots \quad, \frac{d^{n-1} y}{dt^{n-1}}(t_0) = {y^{(n-1)}}_0$$

Again, such an initial value has a unique solution.

Equivalent 1st Order System

An initial value problems for a higher order equation can be reduced to an equivalent initial problem for a system of first order equations. This is important in applications, because there are very few general numerical methods applicable to higher order equations other than those arising from numerical methods for the equivalent first order system. Given an $n$-th order equation

$$\frac{d^ny}{dt^n} = f\left(t, y, \frac{dy}{dt}, \frac{d^2y}{dt^2}, \dots, \frac{d^{n-1}y}{dt^{n-1}} \right)$$ (2)

introduce $n$ new variables

$$w_1 = y$$

$$w_2 = \frac{dy}{dt}$$

$$w_3 = \frac{d^2y}{dt^2}$$

$$\vdots$$

$$w_n = \frac{d^{n-1}y}{dt^{n-1}}$$

There are $n-1$ relationships between the new variables:

$$\frac{dw_1}{dt} = w_2$$

$$\frac{dw_2}{dt} = w_3$$

$$\vdots$$

$$\frac{dw_{n-1}}{dt} = w_n$$

which together with equation (2), which can be rewritten as,

$$\frac{dw_n}{dt} = f(t, w_1, w_2, \dots, w_n)$$

gives a system of $n$ first order equations for the $n$ unknown functions $w_1, \dots, w_n$. The initial values for the $n$-th order equation

$$y(t_0) = y_0, \qquad \frac{dy}{dt}(t_0) = {y'}_0, \quad \dots \quad, \frac{d^{n-1} y}{dt^{n-1}}(t_0) = {y^{(n-1)}}_0$$

are equivalent to the following initial values for our new variables

$$w_1(t_0) = y(t_0) = y_0$$

$$w_2(t_0) = \frac{dy}{dt}(t_0) = {y'}_0$$

$$w_3(t_0) = \frac{d^2y}{dt^2}(t_0) = {y''}_0$$

$$\vdots$$

$$w_n(t_0) = \frac{d^{n-1}y}{dt^{n-1}}(t_0) = y^{(n-1)}_0$$

Thus any initial value problem for an $n$-th order differential equation can be converted to an initial value problem for a system of $n$ first order equations. Similarly, initial value problems for systems of higher order equations can be converted to equivalent problems for 1st order systems using the same techniques. The method is obvious but tedious to write out.

Example:

According to Newton's laws the the position $(x(t), y(t))$ in the orbital plane of a body orbiting a much larger body of mass $M$ centered at the origin is governed by the differential equations

$$\frac{d^2x}{dt^2} = - GM \frac{x}{r^3}$$

$$\frac{d^2y}{dt^2} = - GM \frac{y}{r^3}$$

where $G$ is the gravitational constant and $r = \sqrt{x^2 + y^2}$. Each of the differential equations is 2nd order so the pair of equations is equivalent to a system of four 1st order equations. Introducing variables

$$w_1 = x$$

$$w_2 = \frac{dx}{dt}$$

$$w_3 = y$$

$$w_4 = \frac{dy}{dt}$$

We have the relations

$$\frac{dw_1}{dt} = w_2$$

and

$$\frac{dw_3}{dt} = w_4$$

The original 2nd order equations can be written in terms of the new variables as

$$\frac{dw_2}{dt} = -GM \frac{w_1}{(w_1^2 + w_3^2)^{3/2}}$$

and

$$\frac{dw_4}{dt} = -GM \frac{w_3}{(w_1^2 + w_3^2)^{3/2}}$$

giving our system of four 1st order equations.

Euler's Method

The simplest numerical method for solving initial value problems for differential equations is Euler's method. Suppose we have an initial value problem

$$\frac{dy}{dt} = f(t, y) , \qquad y(t_0) = y_0 .$$

Choose a step-size $\Delta t$, then, by Taylor's approximation

$$y(t_0 + \Delta t) \approx y(t_0) + \frac{dy}{dt}(t_0) \Delta t$$

The differential equation says $dy/dt = f(t, y)$ so that

$$y(t_0 + \Delta t) \approx y_0 + f(t_0, y_0) \Delta t$$

Denote the approximate solution obtained in this way by $\tilde{y}(t)$, then Euler's method can be written

$$\tilde{y}(t_0 + \Delta t) = y_0 + f(t_0, y_0) \Delta t$$

Euler's method has a simple geometric interpretation - we approximate the solution curve through a point $(t_0, y_0)$ by its tangent, which we know has the slope $f(y_0, t_0)$.

We have seen how one step of Euler's method works. The same procedure can be repeated over any number of steps: start with the initial condition

$$\tilde{y}(t_0) = y_0$$

and then apply Euler's method to obtain the approximations

$$\tilde{y}(t_0 + \Delta t) = \tilde{y}(t_0) + f(t_0, \tilde{y}(t_0)) \Delta t$$

$$\tilde{y}(t_0 + 2\Delta t) = \tilde{y}(t_0 + \Delta t) + f(t_0+ \Delta t, \tilde{y}(t_0 + \Delta t)) \Delta t$$

$$\tilde{y}(t_0 + 3\Delta t) = \tilde{y}(t_0 + 2\Delta t) + f(t_0+ 2\Delta t, \tilde{y}(t_0 + 2\Delta t)) \Delta t$$

$$\tilde{y}(t_0 + 4\Delta t) = \tilde{y}(t_0 + 3\Delta t) + f(t_0+ 3\Delta t, \tilde{y}(t_0 + 3\Delta t)) \Delta t$$

$$\vdots$$

Example:

We will use as our example the differential equation

$$\frac{dy}{dt} = -2 t y$$

This equation can be solved by separation of variables to give the general solution

$$y(t) = A e^{-t^2}$$

The initial value problem

$$\frac{dy}{dt} = -2 t y, \qquad y(0) = 1$$

then has the solution

$$y(t) = e^{-t^2}$$

We will write a Scilab function euler(f, t0, y0, dt, n) to approximate the solution of

$$\frac{dt}{dt} = f(t,y)$$

using Euler's method with initial conditions $t$ = t0, $y$ = y0, and time-step dt over n steps.

function [t,y] = euler(f, t0, y0, dt, n)
  t = zeros(n+1,1)
  y = zeros(n+1,1)
  t(1) = t0
  y(1) = y0
  for i = 1:n
    t(i+1) = t(i) + dt
    y(i+1) = y(i) + f(t(i), y(i))*dt
  end
endfunction

Applying euler to our example:

-->function y = func(t, y)
-->  y = - 2*t*y          
-->endfunction  

-->[t,y] = euler(func, 0, 1, 0.01, 200);
 
-->plot2d(t, y)

We can compare the computed solution to the exact solution:

-->err = y - exp(-t.*t);
 
-->plot2d(t, err)

We see that the error is small enough, of the order of $10^{-3}$, for the computed and exact solutions to be barely distinguishable on a graph. The shape of the error curve is not significant. [Aside: the general shape of the error curve can be simply explained. From the graphical interpretation of Euler's method it follows that Euler's method overestimates the solution when the slope of the solution is decreasing (as in the diagram) and underestimates the solution when the slope of the solution is increasing. The change from one type of behaviour to the other will occur near a point of inflection, which, in the exact solution to our example occurs at $x = 1/\sqrt{2}$. There are also other types of error to consider as well to give a complete account.] Euler's method does always perform as well as it does in the example above. Here is an example where it gives wildly inaccurate results:

-->[t,y]=euler(func,0,1,1,10);    
 
-->y'
 ans  =
 
 
         column 1 to 9
 
!   1.    1.  - 1.    3.  - 15.    105.  - 945.    10395.  - 135135. !
 
         column 10 to 11
 
!   2027025.  - 34459425. !

The important change here was the step size. It is clear from the derivation of Euler's method that the accuracy should increase as the step size is decreased, but this example shows that if the step size is too large then totally inaccurate results can be obtained.

Systems of Equations

Euler's method for a single equation takes the form

$$y(t + \Delta t) = y(t) + f(t, y) \Delta t .$$

Simply replacing the scalars $y$ and $f$ by vectors gives a formula which makes sense for systems of equations

$$\mathbf{y}(t + \Delta t) = \mathbf{y}(t) + \mathbf{f}(t, \mathbf{y}) \Delta t .$$

In terms of the individual equations this reads

$$y_1(t + \Delta t) = y_1(t) + f_1(t, y_1, \dots, y_n) \Delta t$$

$$y_2(t + \Delta t) = y_2(t) + f_2(t, y_1, \dots, y_n) \Delta t$$

$$\vdots$$

$$y_n(t + \Delta t) = y_n(t) + f_n(t, y_1, \dots, y_n) \Delta t$$

We will adapt our previous procedure euler to deal with systems of equations. Again it will take the form

     [t, y] = euler(f, t0, y0, dt, n)

and perform n steps of Euler's method with step size dt. The difference now is that since we have a system of equations we need to careful about the sizes of the variables. If we have a system of k equations then we specify:

t:	a column vector of size n+1.
y:	a matrix of size n+1 $\times$ k,
	the rows of y are the components at a given value of $t$,
f(t,y):	the second argument is a row vector of size k,
	it returns a row vector of size k.
y0:	a row vector of size k.

Our function euler is essentially the same as for a single equation, except for the dimensions of the variables involved:

  function [t, y] = euler(f, t0, y0, dt, n)
    k = length(y0)       // this determines the number of equations
    t = zeros(n+1, 1)
    y = zeros(n+1, k)
    t(1) = t0
    y(1,:) = y0
    for i = 1:n
      t(i+1) = t(i) + dt
      y(i+1,:) = y(i,:) + f(t(i), y(i,:))*dt
    end
  endfunction

Example

As an example we will look at the Lotka-Volterra equations which are used to model the dynamics of predator-prey interactions. The prey population is denoted by $y_1$ and the predator population by $y_2$. The differential equations are

$$\frac{dy_1}{dt} = y_1(\alpha_1 - \beta_1 y_2)$$

$$\frac{dy_2}{dt} = y_2(- \alpha_2 + \beta_2 y_1)$$

The parameters $\alpha_1$ and $\alpha_2$ are the natural growth rates of the populations, and the parameters $\beta_1$ and $\beta_2$ determine the interaction between the two populations. For our example we will use the Lotka-Volterra equations with parameter values:

$$\alpha_1 = 1.0, \quad \beta_1 = 0.1$$

$$\alpha_2 = 0.5, \quad \beta_2 = 0.02$$

and initial conditions

$$y_1(0) = 100, \qquad y_2(0) = 10 .$$

Here is the function for the right hand side of our system of equations:

-->a1 = 1;
 
-->a2 = 0.5;
 
-->b1 = 0.1;
 
-->b2 = 0.02;

-->function ydot = lv(t, y)
-->  ydot = zeros(y) 
-->  ydot(1) = y(1)*(a1 - b1*y(2))
-->  ydot(2) = y(2)*(-a2 + b2*y(1))
-->endfunction

By setting the values of the parameters outside the function we could easily change them if we wanted to experiment with different parameter values. By setting ydot = zeros(y) we insure that ydot has the same shape as y, (i.e. if y is a row vector then so is ydot). We will solve the differential equations over the interval $t = [0,25]$ using steps of $\Delta t = 0.01$.

-->[t, y] = euler(lv, 0, [100 10], 0.01, 2500);
 
-->plot2d(t, y)

Note that the plot works because plot2d plots each of the columns of y and we have set things up so that the columns correspond to the two populations $y_1$ and $y_2$. The curve with larger values is $y_1$ the prey population. Both populations show noticeable oscillations but they do not oscillate in phase. This has a biological interpretation. First the prey numbers are decreasing as the predators are increasing - the predators are eating the prey and increasing in number because of the plentiful supply of food. Then the predator numbers decrease as the supply of food continues to decreases. When the predator numbers decrease enough the prey numbers start to increase and the cycle repeats itself. These sorts of oscillations has been observed in natural predator/prey populations. Another useful graph in examples like this is a phase-plane plot when one component of the solution, $y_1$, is plotted against the other, $y_2$, rather than plotting each against $t$.

-->plot2d(y(:,1), y(:,2))

The main thing to be gained from a graph like this is how the two populations vary with one another. There is a clear cyclic pattern to the populations over time. Also apparent is the fact that the maxima and minima of the two populations do not coincide.

About this document ...

AMTH247 Lecture 19

Differential Equations I

This document was generated using the LaTeX2HTML translator Version 2K.1beta (1.61)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split +0 lecture.tex

The translation was initiated by amth247 on 2003-05-27