27-08 — Numerical Methods for ODEs

Phase: 27 — Ordinary Differential Equations | Subject: 27-08 Prerequisites: 27-02 (First-Order ODEs), 27-03 (Second-Order Linear ODEs), basic programming Next subject: 28-01 — Modular Arithmetic Review

Learning Objectives

By the end of this subject, you will be able to:

1. Understand why numerical methods are necessary (most ODEs have no analytical solution)
2. Implement Euler's method and understand its error properties
3. Apply the Runge-Kutta 4 (RK4) method and derive its accuracy
4. Analyse stability of numerical schemes (stiff equations, step size constraints)
5. Choose appropriate step sizes and compare method trade-offs

Core Content

1. Why Numerical Methods?

Most real-world ODEs (weather models, chemical kinetics, orbital mechanics) have no closed-form solution. Numerical methods approximate the solution at discrete points.

Given $y' = f(t, y)$ with $y(t_0) = y_0$, we compute $y(t_0 + h), y(t_0 + 2h), ...$ for step size $h$.

2. Euler's Method

The simplest method: approximate the derivative over a small step.

$$y_{n+1} = y_n + h \cdot f(t_n, y_n)$$

Error: Local truncation error is $O(h^2)$, global error is $O(h)$ — first-order accurate.

Problem: Euler's method is unstable for stiff equations and requires tiny $h$ for reasonable accuracy.

3. Runge-Kutta 4 (RK4)

The workhorse of ODE solving. Uses 4 slope evaluations per step:

$$k_1 = f(t_n, y_n)$$ $$k_2 = f(t_n + h/2, y_n + h k_1/2)$$ $$k_3 = f(t_n + h/2, y_n + h k_2/2)$$ $$k_4 = f(t_n + h, y_n + h k_3)$$

$$y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$

Error: Local truncation error $O(h^5)$, global error $O(h^4)$ — fourth-order accurate.

4. Stability Analysis

A method is stable if errors don't grow uncontrollably. For Euler's method applied to $y' = \lambda y$:

Stability requires $|1 + h\lambda| < 1$ (for $\lambda < 0$).

For RK4: $|1 + h\lambda + (h\lambda)^2/2 + (h\lambda)^3/6 + (h\lambda)^4/24| < 1$.

Stiff equations (some solutions decay very fast, others slowly) require implicit methods (Backward Euler, BDF).

5. Adaptive Step Size

Modern solvers (like scipy.integrate.solve_ivp) adjust $h$ automatically: - Estimate local error using two methods of different order - Increase $h$ when error is small (faster) - Decrease $h$ when error is large (more accurate)

Worked Examples

Example 1: Euler's Method

Problem: Use Euler's method with $h = 0.1$ to approximate $y(0.3)$ for $y' = y$, $y(0) = 1$. (True solution: $y = e^x$, so $y(0.3) = e^{0.3} \approx 1.350$.)

Solution:

$f(t, y) = y$. Euler: $y_{n+1} = y_n + 0.1 \cdot y_n = 1.1 y_n$

• $y_0 = 1$
• $y_1 = 1.1 \cdot 1 = 1.1$
• $y_2 = 1.1 \cdot 1.1 = 1.21$
• $y_3 = 1.1 \cdot 1.21 = 1.331$

Answer: $y(0.3) \approx 1.331$. True value: $e^{0.3} \approx 1.350$. Error: $1.4\%$.

Example 2: RK4 on the Same Problem

Problem: Use RK4 with $h = 0.3$ to approximate $y(0.3)$ for $y' = y$, $y(0) = 1$.

Solution:

Single step from $t_0 = 0$ to $t_1 = 0.3$:

$k_1 = f(0, 1) = 1$ $k_2 = f(0.15, 1 + 0.15 \cdot 1) = f(0.15, 1.15) = 1.15$ $k_3 = f(0.15, 1 + 0.15 \cdot 1.15) = f(0.15, 1.1725) = 1.1725$ $k_4 = f(0.3, 1 + 0.3 \cdot 1.1725) = f(0.3, 1.35175) = 1.35175$

$y_1 = 1 + \frac{0.3}{6}(1 + 2(1.15) + 2(1.1725) + 1.35175)$ $= 1 + 0.05(1 + 2.3 + 2.345 + 1.35175) = 1 + 0.05(6.99675) = 1.3498$

Answer: $y(0.3) \approx 1.3498$. True value: $1.3499$. Error: $0.01\%$ — RK4 with one step beats Euler with 3 steps.

Example 3: Stability Analysis

Problem: Determine the maximum stable step size for Euler's method applied to $y' = -15y$, $y(0) = 1$, for $t \in [0, 1]$.

Solution:

Stability condition: $|1 + h\lambda| < 1$ where $\lambda = -15$.

$|1 - 15h| < 1$ $-1 < 1 - 15h < 1$ $-2 < -15h < 0$ $0 < h < 2/15 \approx 0.133$

Answer: Maximum stable step size for Euler is $h_{max} = 2/15 \approx 0.133$. For $h = 0.1$, we need at least 10 steps. This is a stiff problem — RK4 extends stability to $h_{max} \approx 0.27$.

Practice Problems

Problem 1: Apply Euler's method with $h = 0.2$ to approximate $y(0.4)$ for $y' = x + y$, $y(0) = 1$. Compare with the true solution $y = 2e^x - x - 1$.

Problem 2: For $y' = -2y$ with $y(0) = 1$, show that Euler's method with $h = 1.1$ produces an oscillating unstable solution, while $h = 0.5$ is stable.

Problem 3: Apply one step of RK4 to $y' = x^2 + y^2$, $y(0) = 0$, with $h = 0.5$.

Problem 4: Compare the global error of Euler vs RK4 for $y' = \sin(x)$ with $y(0) = 0$, using $h = 0.1$ over $[0, 1]$.

Problem 5: Explain why stiff ODEs require implicit methods (Backward Euler) rather than explicit methods (Forward Euler, RK4).

Summary

1. Numerical methods approximate ODE solutions when analytical methods fail
2. Euler's method is simple but only first-order accurate ($O(h)$ global error)
3. RK4 is fourth-order accurate ($O(h^4)$) and the standard workhorse
4. Stability limits step size for explicit methods, especially for stiff equations
5. Adaptive step size algorithms balance accuracy and efficiency automatically

Pitfalls

• Choosing a step size based on intuition rather than stability or accuracy requirements: For Euler's method on $y' = -50y$, a step size of $h = 0.05$ "feels" small but violates the stability condition $|1 - 50h| \leq 1$, causing the numerical solution to oscillate and blow up. Always check the stability criterion for the method and the stiffness of the problem before choosing $h$. For stiff ODEs, prefer implicit methods.
• Confusing local truncation error with global error: Local truncation error is the error in a single step assuming exact starting data — it's $\mathcal{O}(h^{p+1})$ for a $p$-th order method. Global error is the accumulated error over the full integration interval — it's $\mathcal{O}(h^p)$. Halving $h$ in Euler (order 1) halves global error, but in RK4 (order 4) it reduces global error by factor 16. Students often cite the local order as if it were global.
• Implementing RK4 with incorrect weighting of the $k_i$ coefficients: The RK4 update is $y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$. Common errors: using $y_n + \frac{h}{4}(k_1+k_2+k_3+k_4)$ (equal weights), forgetting to multiply $k_2$ and $k_3$ by 2, or using the wrong denominator. The specific weights $1/6, 2/6, 2/6, 1/6$ come from Simpson's rule and are essential for the $\mathcal{O}(h^4)$ accuracy.
• Treating all ODEs as non-stiff and using explicit methods blindly: The ODE $y' = -1000y + 1000\sin t$ looks innocent but is stiff — explicit methods need $h < 0.002$ for stability. Using a default explicit solver with moderate $h$ produces garbage. If the solution has components evolving on very different time scales, or if the Jacobian has eigenvalues with large negative real parts, consider implicit methods (backward Euler, trapezoidal, BDF).
• Comparing numerical methods solely by their order without accounting for cost per step: RK4 (4th order, 4 evaluations/step) is typically more efficient than Euler (1st order, 1 eval/step), but comparing RK4 to a 2nd-order method with 2 evaluations per step depends on the accuracy required. For the same number of function evaluations, RK4 often wins — but for problems where only low accuracy is needed, a lower-order method with adaptive step size may be more cost-effective. Always think in terms of work vs. accuracy trade-offs.

Key Terms

• Adaptive step size
• Euler's method
• Explicit method
• Implicit method
• Local truncation error
• RK4 (Runge-Kutta 4)
• Stability analysis
• Stiff equation
• Step size $h$
• Truncation error

Quiz

Q1: The global error of Euler's method is:

A) $O(h^2)$ B) $O(h)$ C) $O(1/h)$ D) $O(h^4)$

Correct: B)

• If you chose B: Correct. Euler's method has local truncation error $O(h^2)$ but global error $O(h)$ because errors accumulate over $O(1/h)$ steps.
• If you chose A: That's the local truncation error. The global error is one order lower.
• If you chose C: Error decreases with smaller $h$, not larger.
• If you chose D: That's RK4's global error order. Euler is much less accurate.

Q2: RK4 requires how many function evaluations per step?

A) 1 B) 2 C) 4 D) 6

Correct: C)

• If you chose C: Correct. RK4 computes $k_1, k_2, k_3, k_4$ — four slope evaluations per step.
• If you chose A: That's Euler's method (one evaluation).
• If you chose B: That's midpoint method (two evaluations).
• If you chose D: Higher-order methods like RK5 or RK6 use more, but RK4 uses exactly 4.

Q3: A stiff ODE is characterised by:

A) Having no solution B) Having solutions with vastly different decay rates, requiring very small step sizes for stability C) Being nonlinear D) Having periodic solutions

Correct: B)

• If you chose B: Correct. Stiff ODEs have some components that decay very rapidly, forcing explicit methods to use tiny $h$ for stability even when the slow component doesn't need it.
• If you chose A: Stiff ODEs do have solutions — they're just hard to compute numerically with explicit methods.
• If you chose C: Stiffness can occur in linear ODEs too (e.g., $y' = \lambda y$ with large negative $\lambda$).
• If you chose D: Periodicity is unrelated to stiffness.

Q4: For Euler's method applied to $y' = \lambda y$, stability requires:

A) $h > 2/|\lambda|$ B) $|1 + h\lambda| < 1$ (for $\lambda < 0$) C) $h < |\lambda|$ D) No constraint — Euler is always stable

Correct: B)

• If you chose B: Correct. The stability region for Euler is the disk $|1 + z| < 1$ in the complex plane where $z = h\lambda$.
• If you chose A: The correct bound is $h < 2/|\lambda|$ for stability, not $h > 2/|\lambda|$.
• If you chose C: The constraint involves both $h$ and $\lambda$ through the product $h\lambda$.
• If you chose D: Euler is conditionally stable. Large $h$ causes exponential error growth for $\lambda < 0$.

Q5: Adaptive step size methods work by:

A) Always using the same $h$ B) Estimating local error and adjusting $h$ to keep error within tolerance C) Using $h = 0$ to avoid errors D) Only solving the ODE once

Correct: B)

• If you chose B: Correct. Adaptive methods compare solutions from two methods of different order, estimate error, and adjust $h$ accordingly.
• If you chose A: That's a fixed-step method. Adaptive methods vary $h$.
• If you chose C: $h = 0$ means no progress. The method adjusts $h$ to be small but non-zero.
• If you chose D: One solve gives one point. Adaptive methods make multiple estimates per step.

Q6: The local truncation error of RK4 is:

A) $O(h)$ B) $O(h^2)$ C) $O(h^4)$ D) $O(h^5)$

Correct: D)

• If you chose D: Correct. RK4 has local truncation error $O(h^5)$ and global error $O(h^4)$.
• If you chose C: That's the global error. Local error is one order higher.
• If you chose A: That's Euler's global error. RK4 is much more accurate.
• If you chose B: That's the local error of the midpoint method.

Q7: For the ODE $y' = -100y + 100\sin(t)$, $y(0) = 0$, Euler with $h = 0.1$ will:

A) Produce an accurate solution B) Be unstable because $h\lambda = -10$ is outside the stability region C) Converge slowly but correctly D) Require no special treatment

Correct: B)

• If you chose B: Correct. $h\lambda = 0.1 \cdot (-100) = -10$, and $|1 - 10| = 9 > 1$, outside Euler's stability region.
• If you chose A: With $h = 0.1$, Euler is unstable for this problem. A smaller $h$ or implicit method is needed.
• If you chose C: The solution oscillates and diverges rather than converging slowly.
• If you chose D: Stiff problems always require special treatment with explicit methods.

Q8: RK4 is preferred over Euler's method because:

A) RK4 always uses fewer steps B) RK4 is fourth-order accurate while Euler is first-order — RK4 achieves much higher accuracy per step C) RK4 never fails for stiff equations D) RK4 requires less computation per step

Correct: B)

• If you chose B: Correct. RK4's fourth-order accuracy means global error decreases as $h^4$ vs Euler's $h$, giving dramatically better accuracy for the same $h$.
• If you chose A: RK4 uses 4 function evaluations per step vs Euler's 1, so it uses more computation per step. But it needs fewer steps for the same accuracy.
• If you chose C: RK4 is still explicit and can be unstable for very stiff equations. Implicit methods are needed for those.
• If you chose D: RK4 uses 4 evaluations per step vs Euler's 1. RK4 costs more per step but needs far fewer steps.

Next Steps

You have completed Phase 27 (Ordinary Differential Equations). Moving to 28-01 — Modular Arithmetic Review to begin Phase 28 (Cryptography Mathematics) with the number-theoretic foundations of modern encryption.