Module 5: Solving Ordinary Differential Equations (Weeks 9-11)
This module focuses on Initial Value Problems (IVPs) for Ordinary Differential Equations
(ODEs). These are problems where we know the rate of change of a system and its
starting state, and we want to predict its future behavior.
Topic 9: Initial Value Problems (IVPs)
The Problem: Given a differential equation dy/dt = f(t, y) and an initial condition y(t₀) =
y₀, find the function y(t) for t > t₀.
• dy/dt = f(t, y) describes the system's dynamics (e.g., a population's growth rate,
a particle's acceleration).
• y(t₀) = y₀ is the starting point.
Why Numerical Methods? Most ODEs encountered in science and engineering cannot
be solved analytically. Numerical methods provide a way to construct a sequence of
points (t₀, y₀), (t₁, y₁), (t₂, y₂), ... that approximates the true solution curve y(t).
Key Notation:
• t_n: The current time point.
• y_n: The approximate value of the solution at t_n.
• h: The step size, the time interval between points (h = t_{n+1} - t_n).
1. Euler's Method
• Core Idea: The simplest possible approach. Assume that the slope f(t, y) is
constant over the small interval [t_n, t_{n+1}] and use it to take a linear step
forward. It's equivalent to using the first term of the Taylor series expansion.
• Formula: y_{n+1} = y_n + h * f(t_n, y_n)
• Pros:
o Extremely simple to understand and implement.
• Cons / Stability Issues:
o Low Accuracy: It is a first-order method. Its local truncation error is
proportional to h², and its global error is proportional to h (written as
O(h)). To double the accuracy, you must halve the step size, which
requires twice as many steps.
o Poor Stability: If the step size h is too large, the method can overshoot the
true solution. For stable systems (where the true solution settles down),
Euler's method can produce a numerical solution that oscillates wildly or
"blows up" (grows without bound). This is a major limitation.
2. Improved Methods (Second-Order)
The weakness of Euler's method is that it only uses the slope information at the
beginning of the interval. Improved methods try to get a better estimate of the average
slope over the interval.
A. Heun's Method (Improved Euler)
• Core Idea: A predictor-corrector method.
1. Predict: Use Euler's method to take a rough, "tentative" step to the end of
the interval.
2. Correct: Average the slope at the beginning of the interval with the slope
at the predicted end point. Use this average slope to take the final, more
accurate step.
• Algorithm:
1. Predictor (Euler step): y_p = y_n + h * f(t_n, y_n)
2. Corrector: y_{n+1} = y_n + (h/2) * [f(t_n, y_n) + f(t_{n+1}, y_p)]
• Analysis: This is a second-order method (O(h²) global error). It is significantly
more accurate and stable than Euler's method for the same step size.
B. The Midpoint Method
• Core Idea: Instead of averaging the slopes at the ends, evaluate the slope at the
midpoint of the interval, which is often a better representation of the average
slope.
• Algorithm:
1. Use Euler's method to find the value at the midpoint: y_mid = y_n + (h/2) *
f(t_n, y_n)
2. Evaluate the slope at this midpoint: k_mid = f(t_n + h/2, y_mid)
3. Use this midpoint slope to take the full step from y_n: y_{n+1} = y_n + h *
k_mid
• Analysis: Also a second-order method (O(h²) global error), comparable in
accuracy and stability to Heun's method.
3. Runge-Kutta Methods (The Workhorse)
• Core Idea: If evaluating the slope at one or two points is good, evaluating it at
several carefully chosen points within the interval must be even better. Runge-
Kutta methods combine multiple slope estimates to achieve very high accuracy.
The Classic 4th-Order Runge-Kutta (RK4)
This is the most famous and widely used method, offering an excellent balance of
accuracy, stability, and computational cost.
• Core Idea: Combine four slope estimates (k1, k2, k3, k4) in a weighted average
to get an exceptionally accurate approximation of the average slope over the
interval.
• Algorithm (The Four "k"s):
1. k1 = f(t_n, y_n) (Slope at the beginning of the interval)
2. k2 = f(t_n + h/2, y_n + (h/2)k1) (Slope at the midpoint, using k1 to get
there)
3. k3 = f(t_n + h/2, y_n + (h/2)k2) (Slope at the midpoint, using the more
accurate k2 to get there)
4. k4 = f(t_n + h, y_n + hk3) (Slope at the end of the interval, using k3 to get
there)
• Final Step (Weighted Average): y_{n+1} = y_n + (h/6) * (k1 + 2k2 + 2k3 + k4)
• Analysis: RK4 is a fourth-order method (O(h⁴) global error). This is incredibly
powerful. Halving the step size h reduces the global error by a factor of 16! This is
why it is the "workhorse" of ODE solvers.
4. Systems of ODEs and Adaptive Step-Size Methods
A. Systems of ODEs
• The Problem: Many real-world systems involve multiple interacting variables,
leading to a system of coupled ODEs: dy₁/dt = f₁(t, y₁, y₂, ..., y_m) dy₂/dt = f₂(t,
y₁, y₂, ..., y_m) ... dy_m/dt = f_m(t, y₁, y₂, ..., y_m)
• The Solution: The methods above generalize directly. We simply treat the
solution y and the functions f as vectors. The formulas for Euler, Heun, and RK4
remain the same, but the additions and multiplications become vector
operations. For example, in RK4, each k becomes a vector of slopes [k1_1, k1_2,
...].
B. Adaptive Step-Size Methods
• The Problem: A fixed step size h is inefficient. When the solution is smooth and
changing slowly, a large h is fine. When the solution is changing rapidly (e.g., a
sudden spike), a small h is required for accuracy.
• The Idea: Let the algorithm automatically adjust the step size h as it proceeds.
• How it Works (Simplified):
1. Take one step of size h with a given method (e.g., RK4) to get y_{n+1}.
2. Take two steps of size h/2 to get another, more accurate estimate,
y*_{n+1}.
3. The difference |y_{n+1} - y*_{n+1}| is an estimate of the local error.
4. If the error is larger than a desired tolerance, decrease h and repeat the
step.
5. If the error is much smaller than the tolerance, increase h for the next
step to save computation time.
• Benefit: Achieves a specified accuracy with the minimum number of function
evaluations, making it highly efficient. Famous implementations include the
RK45 (Runge-Kutta-Fehlberg) method, which uses a 4th-order and a 5th-order
method pair to estimate the error on the fly.