AMTH142 Lecture 21

Differential Equations -- Systems and Higher Order Equations

Contents

• First Order Systems
  • Mathematics
  • Lotka-Volterra Equations
  • Euler's Method
  • Example
  • Scilab
  
  
• Higher Order Equations
  • Formulation
  • Equivalent 1st Order System
  • Example
  
  
• About this document ...

First Order Systems

Mathematics

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 .$$

Lotka-Volterra Equations

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.

Euler's Method

Recall that 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 procedure euler from the previous lecture 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 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.

Scilab

The function ode works the same way for systems of equations as for a single equation

    y = ode(y0, t0, t, f)

Here y0 is assumed to be a column vector and the result, y, contains the components of the solution as rows. This is the opposite of our function euler and we need to take transposes when plotting. Here is our example using ode:

-->t = 0:0.01:25;
 
-->y = ode([100 10]', 0, t, lv);
 
-->plot2d(t', y')                                

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

The results are very similar to those obtained with euler, although the phase plane plot shows an almost exactly closed curve indicating exact periodic behaviour in the solutions. (We would expect the more sophisticated method ode to give more accurate solutions than euler.)

Higher Order Equations

Formulation

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}$-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 convenient since we understand how to deal with the latter numerically. 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)$$ (1)

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 (1), 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. To solve these equations numerically, we will let y denote the vector $(w_1, w_2, w_3, w_4)$. To simplify the computation we can choose units such that $GM = 1$, so the right hand side of system is given by the Scilab function

-->function ydot = orbit(t, y)
-->  ydot = zeros(y)
-->  ydot(1) = y(2)                        
-->  ydot(3) = y(4)                                    
-->  ydot(2) = - y(1)*(y(1)^2 + y(3)^2)^(-3/2)
-->  ydot(4) = - y(3)*(y(1)^2 + y(3)^2)^(-3/2)
-->endfunction

With the initial conditions

$$x(0) = 0.4 \quad \frac{dx}{dt}(0) = 0$$

$$y(0) = 0 \quad \frac{dy}{dt}(0) = 2$$

and solving over the interval $t = [0,50]$:

-->y0 = [0.4 0 0 2]';
 
-->t = 0:0.01:50;
 
-->y = ode(y0, 0, t, orbit);
 
-->plot2d(y(1,:),y(3,:))                      

// position of central body
 
-->plot2d(0,0,-2)

Note that y(1,:) and y(3,:) are the data corresponding to the variables $x(t)$ and $y(t)$. The graph shows that the orbit is an ellipse with the central body at a focus, something well known to Newton.

About this document ...

AMTH142 Lecture 21

Differential Equations -- Systems and Higher Order Equations

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 Applied Mathematics 142 on 2002-10-23