# Predator-Prey Problem

## Contents

and put it in the same directory as your other m-files.

## System of differential equations

Let y1 denote the number of rabbits (prey), let y2 denote the number of foxes (predator). The population growth is described by

• y1' = a1*y1
• y2' = a2*y2

Here the growth rate a1 for rabbits is positive for y2=0, but decreases with increasing y2; the growth rate a2 for foxes is negative for y1=0, but increases with increasing y1:

• a1 = 2 - .5*y2
• a2 = -1 + .5*y1

Hence we obtain the following system of ODEs

• y1' = (2-.5*y2)*y1
• y2' = (-1+.5*y1)*y2

## Given initial conditions

The initial rabbit population is 6, the initial fox population is 2:

• y1(0) = 6
• y2(0) = 2

## Definition of the function f

We have an ODE system y' = f(t,y) where y is a vector and the vector-valued function f is given by

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

## Solve the ODE using the Euler method

Assume that at time t we have the vector y. Then we compute the slope vector s=f(t,y) and define the new y-vector by y = y + h*s.

We repeat this until we reach the final time T.

Note that the values of the Euler method are not very accurate for h=0.1.

```y = [6;2]; t = 0;           % given initial value
h = .1;                     % step size

for i=1:100
plot(t,y(1),'b.',t,y(2),'r.'); hold on    % plot rabbits in blue, foxes in red
s = f(t,y);               % find slope vector
y = y + h*s;              % find new y-vector by using s
t = t + h;
end
hold off
legend('rabbits','foxes')
```

## Solve the initial value problem using ode45

We want to find the functions y1(t) and y2(t) for t in [0,15].

In Matlab we use [ts,ys] = ode45(f,[t0,T],y0).

Here the input arguments are:

• f is the ODE function (can be defined by an @-function, or an m-file)
• t0 is the initial time
• T is the final time
• y0 is the initial y-vector at time t0.

This returns the arrays ts, ys:

• ts is the vector of time values (Matlab chooses h automatically at each step)
• ys is an array containing the values of y: row j contains the values of y1 (rabbits) and y2 (foxes) at the time ts(j) .

[ts,ys] gives a table with t in column 1, y1 (rabbits) in column 2, y2 (foxes) in column 3

```[ts,ys] = ode45(f,[0,15],[6;2]);
[ts,ys]                          % table with values of t, y1, y2 in each row
```
```ans =
0    6.0000    2.0000
0.0251    6.1486    2.1050
0.0502    6.2923    2.2196
0.0754    6.4296    2.3446
0.1005    6.5591    2.4808
0.2069    6.9969    3.1941
0.3133    7.1111    4.1963
0.4197    6.7843    5.4927
0.5262    6.0208    6.9653
0.6137    5.1515    8.1505
0.7013    4.1971    9.1635
0.7888    3.2782    9.9004
0.8763    2.4945   10.3010
0.9560    1.9375   10.3842
1.0357    1.5046   10.2662
1.1153    1.1755    9.9993
1.1950    0.9308    9.6282
1.3057    0.6956    9.0121
1.4163    0.5356    8.3476
1.5270    0.4187    7.6803
1.6377    0.3389    7.0278
1.7386    0.2951    6.4553
1.8395    0.2644    5.9184
1.9404    0.2430    5.4190
2.0413    0.2288    4.9575
2.1513    0.2199    4.4964
2.2612    0.2165    4.0768
2.3711    0.2179    3.6961
2.4811    0.2237    3.3517
2.6374    0.2395    2.9189
2.7937    0.2645    2.5460
2.9499    0.3002    2.2261
3.1062    0.3487    1.9527
3.3371    0.4508    1.6229
3.5679    0.6020    1.3683
3.7988    0.8248    1.1784
4.0296    1.1519    1.0472
4.1967    1.4775    0.9889
4.3638    1.9031    0.9627
4.5309    2.4532    0.9754
4.6980    3.1507    1.0418
4.8650    4.0046    1.1925
5.0321    5.0184    1.4699
5.1992    6.0903    1.9739
5.3663    6.9628    2.8857
5.4157    7.1245    3.2691
5.4652    7.2152    3.7155
5.5146    7.2210    4.2281
5.5640    7.1303    4.8061
5.6135    6.9355    5.4442
5.6629    6.6364    6.1291
5.7124    6.2414    6.8408
5.7618    5.7669    7.5537
5.8406    4.9091    8.6141
5.9193    4.0192    9.4919
5.9981    3.1857   10.1190
6.0769    2.4763   10.4609
6.1518    1.9408   10.5379
6.2268    1.5209   10.4293
6.3017    1.1983   10.1813
6.3767    0.9553    9.8339
6.4820    0.7155    9.2451
6.5873    0.5508    8.6030
6.6926    0.4307    7.9520
6.7980    0.3479    7.3110
6.8961    0.3000    6.7334
6.9942    0.2660    6.1891
7.0924    0.2419    5.6807
7.1905    0.2252    5.2090
7.2970    0.2139    4.7376
7.4035    0.2081    4.3069
7.5101    0.2069    3.9147
7.6166    0.2098    3.5584
7.7638    0.2204    3.1202
7.9110    0.2385    2.7388
8.0582    0.2650    2.4080
8.2054    0.3012    2.1220
8.4423    0.3854    1.7435
8.6792    0.5125    1.4502
8.9161    0.7027    1.2285
9.1530    0.9859    1.0699
9.3284    1.2790    0.9914
9.5038    1.6692    0.9458
9.6791    2.1839    0.9381
9.8545    2.8522    0.9805
9.9937    3.5072    1.0657
10.1329    4.2840    1.2148
10.2721    5.1624    1.4650
10.4112    6.0746    1.8844
10.5057    6.6519    2.3133
10.6002    7.1044    2.9173
10.6947    7.3348    3.7426
10.7892    7.2474    4.8119
10.8574    6.9411    5.7304
10.9256    6.4349    6.7267
10.9938    5.7649    7.7362
11.0620    4.9923    8.6823
11.1302    4.1965    9.4829
11.1984    3.4424   10.0883
11.2666    2.7736   10.4783
11.3348    2.2138   10.6557
11.4057    1.7487   10.6446
11.4766    1.3848   10.4801
11.5475    1.1036   10.2026
11.6184    0.8899    9.8454
11.7113    0.6877    9.3046
11.8041    0.5445    8.7247
11.8970    0.4397    8.1366
11.9899    0.3646    7.5568
12.0920    0.3089    6.9412
12.1941    0.2699    6.3603
12.2962    0.2425    5.8184
12.3983    0.2237    5.3165
12.5019    0.2116    4.8476
12.6055    0.2048    4.4178
12.7090    0.2025    4.0252
12.8126    0.2041    3.6675
12.9561    0.2125    3.2250
13.0996    0.2278    2.8383
13.2430    0.2507    2.5014
13.3865    0.2823    2.2088
13.6312    0.3609    1.7983
13.8759    0.4819    1.4818
14.1207    0.6659    1.2432
14.3654    0.9441    1.0725
14.5240    1.1946    0.9959
14.6827    1.5195    0.9458
14.8413    1.9384    0.9247
15.0000    2.4731    0.9392
```

## Plot both y1(t) and y2(t) vs. t

We can plot the rabbit population vs. t using plot(ts,ys(:,1)) since ys(:,1) gives the first column of the array ys with the rabbit values.

We can plot both rabbits and foxes vs. t using plot(ts,ys). Matlab now plots the first column of ys (rabbits) in blue, and the second column of ys (foxes) in red.

It appears that the functions y1(t), y2(t) are periodic with a period of about 5.

```plot(ts,ys)              % plot rabbit and fox population vs. t
legend('rabbits','foxes'); xlabel('time t')
```

## Plot y2 vs. y1 (phase plane plot)

We can also consider the solution (y1(t),y2(t)) as tracing a curve in the (y1,y2) plane ("phase plane"). We can plot this curve with plot(ys(:,1),ys(:,2)).

The phase plane plot should be a closed curve since the solution is periodic.

The vector field given by f shows the velocity vectors with which the point (y1(t),y2(t)) moves along the trajectory. We see that the point moves along the closed curve counterclockwise as t increases.

Note that we do not quite get a closed curve: after one period we arrive at a slightly different point. This is due to the error of the numerical method in ode45.

```plot(ys(:,1),ys(:,2)); hold on               % plot solution in phase plane
xlabel('rabbits'); ylabel('foxes')
vectfield(f,0:.8:7.5,0:.8:10.5); hold off    % plot vector field
```

## Solve ODE with higher accuracy, then make phase plane plot

We can define additional options using odeset. Here we specify a tolerance of 1e-10 for both the relative error ('RelTol') and the absolute error ('AbsTol').

Now the errors are much smaller, and we obtain a closed curve (at least for the visual accuracy of the plot).

```opt = odeset('RelTol',1e-10,'AbsTol',1e-10);  % use smaller tolerances
[ts,ys] = ode45(f,[0,15],[6;2],opt);          % use opt defined by odeset
plot(ys(:,1),ys(:,2)); hold on                % plot solution in phase plane
xlabel('rabbits'); ylabel('foxes')
vectfield(f,0:.8:7.5,0:.8:10.5); hold off     % plot vector field
```