Section 8.4: Ordinary Differential equations

Created by Tomas de-Camino-Beck
tomasd@math.ualberta.ca

SIR model

Let i(t) be the infected at time t and s(t) suceptibles.  Note that we are not really interested in recovery:

eq1 = s'[
            t] == -β s[t] i[
              t];
eq2 = i'[
            t] == β s[t] i[t] - α i[t];

To find a numreical solution define α and β, then we use NDSolve:

α = 0.04;
β = 0.0002;
solution = 
      NDSolve[{eq1, eq2, s[0] == 990.0, i[0] == 10}, {i, s}, {t, 0, 
          100}]

Plot the numreical solution:

ParametricPlot[
        Evaluate[{t, s[t]} /. solution], {t, 0, 100}];

ParametricPlot[
        Evaluate[{t, i[t]} /. solution], {t, 0, 100}];

This is a plot of s(t) against i(t):

ParametricPlot[
        Evaluate[{i[t], s[t]} /. solution], {t, 0, 100}];

Logistic Equation

Find a solution for the logistic equation

sol = 
      DSolve[{x'[t] == r x[t] (1 - x[t]), x[0] == x0}, 
        x, t]

We can plot a vector field to see the trayectory of solutions:

<< Graphics`PlotField`

r = 2;

note that we ommit the fact that x is a function of time i.e. we use the non-dimensionalized logistic equation:

vectorplot = 
      PlotVectorField[{1, r x (1 - x)}, {t, 0, 6}, {x, 0, 1.3}, 
        AspectRatio -> (0.8/1), 
        Frame -> True, 
FrameTicks -> 
          Automatic]

Plot some solutions

solplot = 
      Plot[Evaluate[{
x[t] /. sol /. 
              x0 -> 0.01, 
x[t] /. sol /. 
              x0 -> 0.2}], {t, 0, 6}, PlotRange -> All, 
        Ticks -> None]

Put them toghether:

Show[vectorplot, solplot]

Predator Prey Model

Consider the following non-dimensional predator v and  prey u model:

prey = u'[t] == 
          u[t] (1 - u[t]/k) - u[t] v[t];
pred = v'[t] == g (u[t] - 1) v[t];

Solving for steady-states, we use the Solve function, letting u' (t), v' (t) = 0:

solution = 
      Solve[{prey /. {u'[t] -> 0}, 
          pred /. {v'[t] -> 0}}, {u[t], v[t]}]

The following code, calculates the Jacobian of a system of equations:

JacobianMatrix[f_List, var_List] := Outer[D, f, var];

Now lets compute the jacobian of the predator prey-system (note that to get the right hanbd side of the equation we use f[[2]]:

J = 
      JacobianMatrix[{prey[[2]], 
          pred[[2]]}, {u[t], v[t]}]; 
    MatrixForm[J]

Evaluated at the steady states:

J /. solution

Note that we get our three matrices, that correspond to the steady states. Now, if we choose some parameter values for g and k, we can compute the eigenvalues in order to study the stability of this system:

param = {g -> 1, k -> 2};

evalJ = J /. solution /. param

In matrix from:

TableForm[MatrixForm /@ evalJ]

We can also study the vector field:

<< Graphics`PlotField`

vp = 
      PlotVectorField[{u (1 - u/k) - u v, 
            g (u - 1) v} /. {g -> 1, k -> 2}, {u, 0, 2}, {v, 
          0, 1}, Frame -> True,  PlotPoints -> 15]

solution1 = 
      NDSolve[{prey /. param, pred /. param, u[0] == 2, 
          v[0] == 0.1}, {u, v}, {t, 0, 100}]

phase1 = 
        ParametricPlot[
          Evaluate[{u[t], v[t]} /. solution1], {t, 0, 100}];

Show[vp, phase1]