Systems of differential equations

Ivan Cnop
Vrije universiteit Brussel
icnop@vub.ac.be

Examples of systems of differential equations

This notebook offers templates for solving linear systems of differential equations.
The examples are classified by the entries of the coefficient matrix.
Students should repeat the computations for many initial values and slightly different coefficient sets.

Communicating containers with feedback

Coefficient matrix

a = -.4; b = .1; c = .4; d = -.2;

A=(({{a, b}, {c, d}}))

A // MatrixForm

Eigensystem[A]

Solution

system = {x'[t] == First[A] . {x[t], y[t]}, y'[t] == 
          Last[A] . {x[t], y[t]}}

initialValues = {x[0] == 10, y[0] == 0}

sol = 
      DSolve[Flatten[{system, initialValues}], {x[t], y[t]}, t]

Except for the number   e   in the solution, we cannot guess what the solution looks like.

Plane plot

ParametricPlot[{x[t] /. First[First[sol]], y[t] /. Last[First[sol]]}, {
        t, 0, 50}, PlotRange -> All, AxesLabel -> {vat1, vat2}]

Plot[(y[t] /. Last[First[sol]]), {t, 0, 50}, PlotRange -> All]

Competing species

Coefficient matrix

a = .1; b = -.1; c = -.2; d = .4;

A=(({{a, b}, {c, d}}))

A // MatrixForm

Eigensystem[A]

system = {x'[t] == First[A] . {x[t], y[t]}, y'[t] == 
          Last[A] . {x[t], y[t]}}

initialValues = {x[0] == 150, y[0] == 200}

Solution

sol = 
        DSolve[Flatten[{system, initialValues}], {x[t], y[t]}, 
          t];

Plane plot

ParametricPlot[{x[t] /. First[First[sol]], y[t] /. Last[First[sol]]}, {
        t, 0, 5}, PlotRange -> All, AspectRatio -> Automatic]

ParametricPlot[{x[t] /. First[First[sol]], y[t] /. Last[First[sol]]}, {
        t, 0, 20}, PlotRange -> All]

The population   x   dies out first.

Plot[{x[t] /. First[First[sol]], y[t] /. Last[First[sol]]}, {t, 0, 15}, 
      PlotRange -> {0, 1000}, PlotStyle -> {Hue[0], Hue[.7]}]

Now   x   is going to zero rapidly.

Different initial values

initialValues = {x[0] == 400, y[0] == 200}

sol = 
        DSolve[Flatten[{system, initialValues}], {x[t], y[t]}, 
          t];

Plane plot

ParametricPlot[{x[t] /. First[First[sol]], y[t] /. Last[First[sol]]}, {
        t, 0, 10}, AxesOrigin -> {0, 0}]

ParametricPlot[{x[t] /. First[First[sol]], y[t] /. Last[First[sol]]}, {
        t, 0, 20}, PlotRange -> All]

Again, the population   x   dies out first.

Plot[{x[t] /. First[First[sol]], y[t] /. Last[First[sol]]}, {t, 0, 15}, 
      PlotRange -> {0, 1000}, PlotStyle -> {Hue[0], Hue[.7]}]

Here the population   y   dies out first.

Conclusion:

In this example the initial value decides who wins.

Look for initial values that ensure both  x   and   y   a reasonably long life span.

A mechanical example: pendulum with friction

Coefficient matrix

Friction is proportional to velocity, and acts against it:

a = 0; b = 1; c = -1; d = -.2;

A=(({{a, b}, {c, d}}))

A // MatrixForm

Eigensystem[A]

system = {x'[t] == First[A] . {x[t], y[t]}, y'[t] == 
          Last[A] . {x[t], y[t]}}

initialValues = {x[0] == 10, y[0] == 0}

Solution

sol = 
        DSolve[Flatten[{system, initialValues}], {x[t], y[t]}, 
          t];

Phase plane plot

ParametricPlot[{x[t] /. First[First[sol]], y[t] /. Last[First[sol]]}, {
        t, 0, 50}, PlotRange -> All]

Plot position (red) and velocity (blue) together:

Plot[{x[t] /. First[First[sol]], y[t] /. Last[First[sol]]}, {t, 0, 15}, 
      PlotRange -> {-10, 10}, PlotStyle -> {Hue[0], Hue[.7]}]

A medical application

Coefficient matrix

Here  x  is the glucose released by the liver (above normal values ) and  y   is the extra level of insuline produced by the pancreas that will reduce the glucose level.

a = -.1; b = -.2; c = .3; d = -.2;

A=(({{a, b}, {c, d}}))

A // MatrixForm

Eigensystem[A]

system = {x'[t] == First[A] . {x[t], y[t]}, y'[t] == 
          Last[A] . {x[t], y[t]}}

Solution

initialValues = {x[0] == 1, y[0] == 0}

sol = 
        DSolve[Flatten[{system, initialValues}], {x[t], y[t]}, 
          t];

Plane plot

\(ParametricPlot[

      	{x[t] /. First[First[sol]], y[t] /. Last[First[sol]]}, {t, 0, 50}, 
      PlotRange -> All, AxesLabel -> {glucose, insuline}]\)

This system stabilises. Here is a plot of the extra  extra insuline level in time:

Plot[(y[t] /. Last[First[sol]]), {t, 0, 50}]

Predator-prey system

Coefficient matrix

Here   x   is the large number of rats reproducing rapidly and   y   is the small number of owls feeding on rats.

a = 2.1; b = -0.104; c =  .4; d = 1.5;

A=(({{a, b}, {c, d}}))

A // MatrixForm

Eigensystem[A]

system = {x'[t] == First[A] . {x[t], y[t]}, y'[t] == 
          Last[A] . {x[t], y[t]}}

Solution

initialValues = {x[0] == 10000, y[0] == 10}

sol = 
        DSolve[Flatten[{system, initialValues}], {x[t], y[t]}, 
          t];

Plane plot

ParametricPlot[{x[t] /. First[First[sol]], 
        y[t] /. Last[First[sol]]}, {t, 0, 4}, PlotRange -> All, 
      AxesLabel -> {rats, owls}]

ParametricPlot[{x[t] /. First[First[sol]], 
        y[t] /. Last[First[sol]]}, {t, 0, 100}, PlotRange -> All, 
      AxesLabel -> {rats, owls}]

The ratio tends to a constant:

Plot[(y[t] /. Last[First[sol]])/(x[t] /. 
            First[First[sol]]), {t, 0, 20}, PlotRange -> All]