Description

1. Consider a Runge–Kutta method starting at (𝑡𝑛,𝑥𝑛) that takes an Euler step of length ℎ/2 to (𝑡𝑛+1/2,𝑥𝑛+1/2) and then uses the new evaluation at that point to take a complete Euler step from (𝑡𝑛,𝑥𝑛) of length ℎ.
Find the order of this method and write down its Butcher tableau. We will refer to it as the “midpoint method”.
2. Define a type RKMethod to represent a general explicit Runge–Kutta method defined by a Butcher tableau as follows:
struct RKMethod{T} c::Vector{T} b::Vector{T} a::Matrix{T} s::Int # number of stages end
Make it into a function by completing the function
function (method::RKMethod)(f, x, t, h) … end
to execute one step of the corresponding Runge–Kutta method with initial condition 𝑥 at time 𝑡 and step size ℎ. Your code should work for both scalar and vector 𝑥 and a possibly vector-valued function 𝑓 = 𝑓(𝑡,𝑥)
(Assume that 𝑎 is a lower-triangular matrix, corresponding to an explicit method.)
3. Define RK methods euler, midpoint and RK4 using their respective tableaus.
4. Write a routine integrate with the signature function integrate(method, f, x0, t0, t_final, h) where method is a RK method as defined above and ℎ is a fixed step size.
Make sure that the final step lands exactly at the final time by taking that final step as a special case. 5. Use each method on the ODE 𝑥̇ = 1.5𝑥 with 𝑥0 = 2 and integrate from 𝑡 = 0 for a time 𝑡final = 3.
Find the rate of convergence of the numerical solution to the exact solution as ℎ → 0 for each method. Do they correspond to our analytical expectations?
𝑢̇(𝑡) = exp [𝑡 − 𝑢 sin(𝑢)].
Integrate it using RK4 from 𝑡 = 0 to 𝑡 = 5 with a step size ℎ = 10−2. Now integrate it with a step size ℎ = 10−3.
Plot both solutions 𝑥(𝑡) as a function of 𝑡. What do you observe? What do you think is happening?
Exercise 2: Adaptivity in the Euler method
In this exercise we will invetigate adaptivity in ODE solvers by taking the simplest case: an adaptive Euler method.
1. Consider one step of the Euler method. Write down the local (single-step) error in terms of the step size ℎ and the unknown constant 𝐶. Call the approximation obtained at the end of the step 𝑥(1).
2. Now consider taking two consecutive Euler steps of size ℎ/2. Would you expect this to give a better or a worse approximation to the true solution? Write down the total error after taking the two steps, assuming that the constant 𝐶 is the same for both. Call the approximation at the end of this combined step 𝑥(2).
3. Define Δ𝑥 as the difference between the two approximations. Use this to find the step size ℎ′ that will give an error per unit time of a given size
𝜖.
4. Use this derivation to write an adaptive Euler integrator adaptive_euler(f, t0, t_final, epsilon) that tries to keep the global error less than 𝜖, using an update rule similar to the one we discussed in class. Add a multiplicative factor 0.9 to that rule to make the method behave better.
5. Use it to integrate the same ODE as in exercise 1. Plot the step size taken as a function of time.
6. Now integrate the equation for the van der Pol oscillator:
𝑥̈ − 𝜇(1 − 𝑥2)𝑥̇ + 𝑥 = 0,
with 𝜇 = 5. Use initial conditions 𝑥0 = 0 and 𝑥̇0 = 1.
7. Plot trajectories in (𝑥,𝑥)̇ phase space and (separately) the solution 𝑥(𝑡) as a function of 𝑡.
8. Plot the step size as a function of time. What do you observe? How do you interpret this?
Exercise 3: SIR model In this exercise we will study the SIR model of the dynamics of an infectious disease outbreak (“epidemic”) in a population.
1. Use e.g. RK4 to solve the SIR equations:
𝑆 = −𝛽𝑆 𝐼̇ (1)
𝐼 = 𝛽𝑆 𝐼 − 𝛾𝐼̇ (2)
𝑅 = 𝛾𝐼̇ (3)
Here 𝐼 is the proportion of the population which is infectious. 𝛽 is the rate of contact between susceptible and infectious individuals, and 𝛾 the recovery rate.
Use 𝑆0 = 0.99 and 𝐼0 = 0.01.
2. Make an interactive visualization, varying 𝛽 and 𝛾 in, say, the range 0 to 1.
3. What do you observe? Can you interpret this?
4. By summing the equations we see that 𝑆 + 𝐼 + 𝑅 should be constant (equal to the total population, assuming no births or deaths). For representative parameter values 𝛽 = 0.1 and 𝛾 = 0.05, how well does the numerics conserve this quantity?