Systems of IVPs
Name:
Approx. time spent on this assignment:
Did you do the suggested reading? Yes/No
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. 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.
Hint: You will want to define
\begin{align} u(t) &= \begin{pmatrix} x(t) \\ y(t) \end{pmatrix}. \nonumber \end{align}
Answer.
Your answer here
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.
Your answer here
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.
Your answer here
Exercise 4. What happens when (x_0, y_0) = (C/D, A/B)? Why?
Answer.
Your answer here
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.
Your answer here
Exercise 6. Write down the system of differential equations governing the variable
\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.
Your answer here
Exercise 7. Use the midpoint method to implement the numerical solution to this system of equations. Plot your solution and comment on the dependence of the solution on \gamma. Hint: You can use the code samples given above.
# Your answer hereRead: Sections 5.1 - 5.4 of Burden, Faires, and Burden (2015)