Example for system of ODEs

clearvars; clf; format short g

You first have to download vectorfield.m and put this file in the folder where your mlx files are.

Population of rabbits and foxes

Let

denote the population of rabbits, let

denote the population of foxes at time t.

We have for the rates of change

Rabbits (prey)

where

Foxes (predator)

where

We use the vector valued function

with derivative

We can write the system of ODEs as

Here we have

In Matlab we can define this function as

f = @(t,y) [ (2-y(2))*y(1) ; (-1+y(1))*y(2) ];

We can evaluate e.g.

using

s = f(0,[3;4])
s = 2×1
    -6
     8

Initial value problem

At the initial time

we have initial populations

for rabbits and foxes, i.e.

Euler method

We want to approximate

using one step of the Euler method with

h = 1/2;

At the initial time

we have the value

t0 = 0;
y0 = [2;2];

We now evaluate

s = f(t0,y0)
s = 2×1
     0
     2

At the time

we obtain the approximation

y1 = y0 + h*s
y1 = 2×1
     2
     3

and obtain

. This means

Improved Euler method

We now perform one step of the improved Euler method with

s1 = f(t0,y0)
s1 = 2×1
     0
     2

yE = y0+h*s
yE = 2×1
     2
     3

s2 = f(t0+h,yE)
s2 = 2×1
    -2
     3

y1 = y0 + h*(s1+s2)/2
y1 = 2×1
          1.5
         3.25

Using ode45 to solve the initial value problem

We want to approximate the solution

for

with

[ts,ys] = ode45(f,[0,9],[2;2]);
disp([ts,ys]) % shows t, y1(t), y2(t) in each row (scroll to see all rows)
            0            2            2
     0.050238       1.9949       2.1029
      0.10048       1.9793       2.2099
      0.15071       1.9532       2.3199
      0.20095       1.9166       2.4315
       0.3277       1.7821       2.7103
      0.45444       1.6029       2.9598
      0.58119       1.4014        3.154
      0.70794       1.1998       3.2763
      0.82671       1.0276         3.32
      0.94548      0.87877       3.3009
       1.0642      0.75606       3.2293
        1.183      0.65797       3.1177
       1.3018      0.58086       2.9793
       1.4206      0.52181       2.8242
       1.5393      0.47768        2.661
       1.6581      0.44592        2.496
       1.8173       0.4193       2.2799
       1.9765      0.40769       2.0761
       2.1358      0.40878       1.8892
        2.295      0.42162       1.7213
       2.4589        0.447       1.5685
       2.6229      0.48512       1.4368
       2.7868      0.53706       1.3259
       2.9507      0.60443       1.2355
       3.1757      0.72492       1.1453
       3.4007      0.88467       1.0952
       3.6257       1.0871       1.0903
       3.8507       1.3278       1.1414
        4.024       1.5305       1.2287
       4.1972       1.7287       1.3711
       4.3705       1.8944       1.5804
       4.5437       1.9907       1.8628
       4.6566       1.9971       2.0851
       4.7694       1.9515       2.3284
       4.8822       1.8544       2.5786
        4.995       1.7137       2.8179
       5.1078       1.5439       3.0254
       5.2206       1.3621       3.1844
       5.3335       1.1833       3.2849
       5.4463       1.0198       3.3231
       5.6029      0.82971       3.2825
       5.7596      0.68451       3.1581
       5.9163      0.58156       2.9778
        6.073      0.51002       2.7696
       6.2055      0.46618       2.5875
        6.338      0.43653       2.4057
       6.4705      0.41849       2.2297
        6.603       0.4104       2.0632
       6.7581      0.41222       1.8828
       6.9132      0.42521       1.7202
       7.0684      0.44915       1.5762
       7.2235      0.48445       1.4509
       7.4205      0.54703       1.3184
       7.6175      0.63236       1.2155
       7.8145      0.74384       1.1425
       8.0115      0.88476        1.101
        8.191       1.0404       1.0934
       8.3705       1.2223       1.1191
         8.55       1.4248       1.1855
       8.7295       1.6331       1.3034
       8.7971       1.7085       1.3639
       8.8647       1.7795       1.4343
       8.9324       1.8439       1.5153
            9       1.8996       1.6074

Here ts is a vector of t-values, starting with

and ending with T.

ys is an array with two columns: the first column contains the approximations for

, the second column contains the approximations for

The last row of ys contains the approximations for the final time T:

ys(end,:)         % last row of array ys, contains values for time T
ans = 1×2
       1.8996       1.6074

This means that

The first column ys(:,1) contains the population of rabbits. We can plot this using plot(ts,ys(:,1))

plot(ts,ys(:,1))

xlabel('t'); title('population y_1(t) of rabbits')

The second column ys(:,2) contains the population of foxes.

We can plot both functions

and

together using plot(ts,ys)

plot(ts,ys); legend('y_1','y_2'); xlabel('t')

title('y_1(t) (rabbits) and y_2(t) (foxes)')

We start at a maximum of rabbits, then the foxes increase to a maximum, then the rabbits decrease to a minumum, then the foxes decrease to a minimum, then the rabbits increase to a maximum, etc.

Note that the solution is periodic: after

everything repeats

Solution in the "phase plane"

The solution

describes a point moving around in the

plane, the so-called phase plane.

The velocity vector of this moving point at time t is given by the vector

In our example the function

does not depend on t, i.e., we have an autonomous ODE.

We can plot the vector

at a grid of points in the

plane, this gives a "vector field":

vectorfield(f,0:.5:3,0:.5:4.25); hold on

xlabel('y_1 (rabbits)'); ylabel('y_2 (foxes)')

Recall that ys(:,1) contains the rabbit population

and ys(:,2) contains the fox population

We can plot plot the points

in the phase plane using plot(ys(:,1),ys(:,2))

[ts,ys] = ode45(f,[0 9],[2;1]);

plot(ys(:,1),ys(:,2),'b')

hold off

The periodic solution

corresponds to a closed curve in the phase plane.