Systems of IVPs
In classes, we considered the following problem: seek u:[0,T] \to \mathbb R such that
\begin{align} u(0) &= u_0 \nonumber\\ u'(t) &= f\big( t, u(t) \big). \tag{IVP} \end{align}
This equation governs the dynamics of a single unknown as it evolves. In practice, there are often multiple unknowns that interact over time and one wishes to model the dynamics as the system evolves.
For example, one may consider the so-called Lotka–Volterra predator–prey model describing the evolution of a biological system involving two species: a prey population of size x(t) at time t, and a predator population of size y(t) at time t. These populations are modelled by the following system of equations: x(0) = x_0, y(0) = y_0 for some initial data (x_0,y_0) and
\begin{align} x'(t) &= A x(t) - B x(t) y(t) \nonumber\\ y'(t) &= -C y(t) + D x(t) y(t) \tag{LV} \end{align}
for some constants A, B, C, D \geq 0
Exercise 1. Write (\text{LV}) as a problem of the form: seek u : [0,T] \to \mathbb R^2 such that
\begin{align} u(0) &= u_0 \nonumber\\ u'(t) &= f\big(t, u(t)\big) \end{align}
for some f: [0,T] \times \mathbb R \to \mathbb R^2.
Answer.
On defining u_1(t) := x(t) and u_2(t) := y(t), we find that
\begin{align} u(0) &= \begin{pmatrix}x_0\\y_0\end{pmatrix} \nonumber\\ u'(t) &= \begin{pmatrix} A u_1(t) - B u_1(t) u_2(t) \\ -C u_2(t) + D u_1(t) u_2(t) \end{pmatrix}. \end{align}
Notice that you may apply the numerical schemes from lectures to this system of IVPs. For example, the midpoint method reads:
function midpoint( u0, f, T, n )
h = T/n
t = 0:h:T
u = fill(float(u0), n+1)
u[1] = u0
for j = 1:n
u[j+1] = u[j] + h * f( t[j] + h/2, u[j] + (h/2)*f(t[j], u[j]) )
end
return u
end midpoint (generic function with 1 method)Here, is the numerical solution of (\text{VP}) with a particular choice of parameters A, B, C, D, initial condition u_0, and mesh \{ t_j \}:
A, B, C, D = 1., 1., .5, 1.
u0 = [1.,.01]
f(t, u) = [A*u[1]-B*u[1]*u[2], -C*u[2]+D*u[1]*u[2]]
T, n = 50., 500
u = midpoint( u0, f, T, n)
u = [U[j] for U ∈ u, j ∈ 1:2]
plot(0:T/n:T, u,
label=["prey" "predator"], l=(1, :black), m=2)
Exercise 2. Explain this plot with reference to the physical system that is being modelled.
Answer.
Initially, there are very few predators so the prey are free to grow. Eventually, this means there is available food for the predators and so the predator population grows as a result. Then, the large size of the predator population needs a large amount of food and so the prey is eaten and the prey population decreases rapidly. This reduces the available food for the predators and so the predator population also drops. This behaviour seems to converge to a stable periodic solution.
Another way to see the dynamics is to plot (x(t), y(t)) in phase space:
X, Y = [], []
anim = @animate for t ∈ 1:n+1
x,y = u[t,:]
plot(X, Y,
title="LV numerical solution",
xlabel=L"x(t)", ylabel=L"y(t)", legend=false)
push!(X, x)
push!(Y, y)
scatter!([x], [y], m=(5, :red))
end
mp4(anim, "pics/LV.gif")[ Info: Saved animation to c:\Users\math5\Math 5485\Math5486\pics\LV.gif
Exercise 3. Play around with the parameters in this model (A, B, C, D). What do you notice about the numerical solution?
Answer.
There are many possible answers here - you should have played around with the following parameters and noticed that the effect makes sense when compared to the physical meaning of the variables:
Exercise 4. What happens when (x_0, y_0) = (C/D, A/B)? Why?
Answer.
In this case, we have u'(t) = 0 and so u(t) = \begin{pmatrix} C/D \\ A/B \end{pmatrix} is a constant solution.
Consider the follwing initial value problem: seek \theta: [0,1] \to \mathbb R such that
\begin{align} \theta(0) &= \theta_0 \nonumber\\ \theta'(0) &= 0 \nonumber\\ \theta''(t) &= - \gamma \theta'(t) -\frac{g}{L} \sin \theta(t). \nonumber \end{align}
This is a second order IVP modelling the motion of a pendulum of length L and with damping \gamma > 0. Here, \theta is the angle made with the vertical position and g \approx 9.81 m/s^2 is the gravitational constant. See the beginning of Chapter 5 of Burden, Faires, and Burden (2015) for a picture.
Exercise 5. Show that \theta(t) = 0 is a solution when \theta_0 = 0. Explain why this makes sense with reference to the physical system we are modelling.
Answer.
If \theta(t) = 0, we have \theta''(t) = \theta'(\theta) = 0 and so - \gamma \theta'(t) -\frac{g}{L} \sin \theta(t) = 0. When the initial position of the pendulum is vertical, the pendulum will stay at rest.
Exercise 6. Write down the system of differential equations governing the quantity
\begin{align} u(t) = \begin{pmatrix} u_1(t) \\ u_2(t) \end{pmatrix} = \begin{pmatrix} \theta(t) \\ \theta'(t) \end{pmatrix}. \end{align}
Answer.
We have
\begin{align} u'(t) = \begin{pmatrix} \theta'(t) \\ \theta''(t) \end{pmatrix} % = \begin{pmatrix} u_2(t) \\ - \gamma u_2(t) -\frac{g}{L} \sin u_1(t) \end{pmatrix} \end{align}
and u(0) = \begin{pmatrix} \theta_0 \\ 0 \end{pmatrix}.
Exercise 7. Use the midpoint method to implement the numerical solution to this system of equations. Plot your solution.
u0, g, γ, L = [1,0], 9.81, .25, 10.
f2(t, u) = [u[2], -γ*u[2] - (g/L) * sin( u[1] )]
T, n = 50., 500
u = midpoint( u0, f2, T, n)
u = [U[j] for U ∈ u, j ∈ 1:2]
plot(0:T/n:T, u[:,1], legend=false,
xlabel=L"t", ylabel=L"\theta(t)",
l=(1, :black), m=2)
Read: Sections 5.1 - 5.4 of Burden, Faires, and Burden (2015)