In Lesson 1 of this chapter, students saw Euler and improved Euler methods, and learned the meaning of local truncation error and global error. We now generalize that idea to the Runge–Kutta family, which achieves higher accuracy without computing higher derivatives explicitly. This is especially useful in system dynamics, where state models are often available only as first-order right-hand sides \( \dot{\mathbf{x}} = \mathbf{f}(t,\mathbf{x}) \).

Throughout the lesson, we study the initial value problem (IVP)

\[ \dot{y}(t) = f(t,y), \qquad y(t_0)=y_0 \]

and later extend the formulas componentwise to vector systems \( \mathbf{y} \in \mathbb{R}^m \). We assume \( f \) is sufficiently smooth so that Taylor expansions are valid.

An explicit \( s \)-stage Runge–Kutta method computes stage slopes \( k_i \) and combines them linearly:

\[ k_i = f\!\left(t_n + c_i h,\; y_n + h\sum_{j=1}^{i-1} a_{ij} k_j\right), \quad i=1,\dots,s \]

\[ y_{n+1} = y_n + h\sum_{i=1}^{s} b_i k_i \]

The coefficients are typically organized in a Butcher tableau. For explicit methods, \( a_{ij}=0 \) whenever \( j \ge i \), so each stage depends only on previously computed stages. This preserves an easy computational pipeline.

The design objective is to choose \( a_{ij}, b_i, c_i \) so that the one-step map matches the Taylor expansion of the exact solution up to a desired order.

flowchart TD
  A["Known state y_n at time t_n"] --> B["Compute stage k1 from f(t_n, y_n)"]
  B --> C["Compute intermediate stages k2, k3, ..., ks"]
  C --> D["Weighted sum of stages"]
  D --> E["Advance to y_(n+1)"]
  E --> F["Repeat with next step size h"]
        

Consider the midpoint RK2 method:

\[ k_1 = f(t_n,y_n), \qquad k_2 = f\!\left(t_n + \frac{h}{2},\; y_n + \frac{h}{2}k_1\right) \]

\[ y_{n+1} = y_n + h k_2 \]

To show it is second order, expand the exact solution around \( t_n \):

\[ y(t_n+h) = y_n + h y'_n + \frac{h^2}{2} y''_n + O(h^3) \]

Since \( y'_n=f(t_n,y_n) \), we expand \( k_2 \) about \( (t_n,y_n) \):

\[ k_2 = f_n + \frac{h}{2} f_t\big|_n + \frac{h}{2} f_y\big|_n f_n + O(h^2) \]

Multiplying by \( h \) gives

\[ y_{n+1} = y_n + h f_n + \frac{h^2}{2}\left(f_t + f_y f\right)\Big|_n + O(h^3) \]

but \( y'' = f_t + f_y f \) along the exact trajectory, hence the method matches the exact Taylor expansion up to \( O(h^2) \) and has second-order global accuracy.

The classical four-stage RK4 method is

\[ \begin{aligned} k_1 &= f(t_n,y_n) \\ k_2 &= f\!\left(t_n+\frac{h}{2},\; y_n + \frac{h}{2}k_1\right) \\ k_3 &= f\!\left(t_n+\frac{h}{2},\; y_n + \frac{h}{2}k_2\right) \\ k_4 &= f\!\left(t_n+h,\; y_n + h k_3\right) \end{aligned} \]

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

A full proof of fourth-order accuracy is long if expanded directly, so numerical analysis uses order conditions. For order 1 through 4, representative conditions are:

\[ \sum_i b_i = 1, \qquad \sum_i b_i c_i = \frac{1}{2}, \qquad \sum_i b_i c_i^2 = \frac{1}{3} \]

\[ \sum_{i,j} b_i a_{ij} c_j = \frac{1}{6}, \qquad \sum_i b_i c_i^3 = \frac{1}{4}, \qquad \sum_{i,j} b_i c_i a_{ij} c_j = \frac{1}{8} \]

The classical RK4 coefficients satisfy the complete set of fourth-order conditions, implying local truncation error \( O(h^5) \) and global error \( O(h^4) \) for smooth problems.

In practice, RK4 is a strong default fixed-step method for non-stiff system dynamics simulations because it significantly improves accuracy over Euler-type methods at modest extra cost (4 function evaluations per step).

Let \( \Phi_h \) denote the exact one-step flow and \( \Psi_h \) the numerical one-step map. A method has order \( p \) if its local defect satisfies

\[ \left|\Phi_h(y_n) - \Psi_h(y_n)\right| \le C h^{p+1} \]

Assuming \( \Psi_h \) is Lipschitz in \( y \) with constant \( 1 + Lh \), the global error \( e_n = y(t_n)-y_n \) satisfies

\[ |e_{n+1}| \le (1+Lh)|e_n| + C h^{p+1} \]

Iterating the inequality (discrete Gronwall) over \( N=(T-t_0)/h \) steps yields

\[ |e_N| \le \frac{C}{L}\left(e^{L(T-t_0)}-1\right) h^p \]

Therefore, a local error of order \( p+1 \) implies global error of order \( p \). This is why RK4 (local \( O(h^5) \)) is globally fourth order.

For a method of global order \( p \), the asymptotic global error often behaves like

\[ E(h) \approx K h^p \]

for some problem-dependent constant \( K \). If a target accuracy \( E_{\mathrm{tol}} \) is prescribed, a first estimate is

\[ h \approx \left(\frac{E_{\mathrm{tol}}}{K}\right)^{1/p} \]

Since \( K \) is unknown, a common engineering workflow is:

1. Run with a baseline step \( h \).
2. Repeat with \( h/2 \).
3. Estimate observed order by comparing solutions.
4. Decrease \( h \) until the change is below tolerance.

For RK4, halving \( h \) should reduce global error by approximately \( 2^4=16 \) in the asymptotic regime.

Fixed-step methods are inefficient when the solution changes slowly in some regions and rapidly in others. Adaptive RK methods estimate local error and choose the next step size automatically.

In this lesson, we use a rigorous and easy-to-implement approach: RK4 with step doubling. Compute

• one full RK4 step of size \( h \) to get \( y_{\text{full}} \),
• two RK4 half-steps of size \( h/2 \) to get \( y_{\text{half2}} \).

Because RK4 has local error \( O(h^5) \), the difference scales as

\[ y_{\text{half2}} - y_{\text{full}} = \left(1-2^{-4}\right) C h^5 + O(h^6) = \frac{15}{16} C h^5 + O(h^6) \]

Hence a local error estimate for the full-step RK4 solution is

\[ e_{\mathrm{loc}} \approx \frac{|y_{\text{half2}} - y_{\text{full}}|}{15} \]

For vector systems, a normalized componentwise criterion is standard. Define

\[ \mathrm{sc}_i = \mathrm{atol} + \mathrm{rtol}\,\max\!\left(|y_i|,|\hat{y}_i|\right) \]

\[ \mathrm{err} = \sqrt{\frac{1}{m}\sum_{i=1}^{m}\left(\frac{e_i}{\mathrm{sc}_i}\right)^2} \]

Accept the step if \( \mathrm{err} <= 1 \). Then update \( h \) using the controller

\[ h_{\text{new}} = s\, h \left(\frac{1}{\mathrm{err}}\right)^{1/5} \]

where \( s \in (0,1) \) is a safety factor (e.g., 0.8 to 0.95). The exponent \( 1/5 \) comes from the \( h^5 \) local error behavior of RK4.

flowchart TD
  A["Start step with current h"] --> B["Compute y_full with one RK4 step"]
  B --> C["Compute y_half2 with two half RK4 steps"]
  C --> D["Estimate local error from difference"]
  D --> E["Is normalized error <= 1 ?"]
  E -->|Yes| F["Accept step and advance state"]
  E -->|No| G["Reject step"]
  F --> H["Update h using safety controller"]
  G --> H
  H --> I["Next trial step"]
        

In control engineering models, step-size selection is not only a numerical-analysis issue; it also affects interpretation of transients, peaks, and switching-like behavior (formal event handling appears in Lesson 4). Practical guidance:

• Choose \( \mathrm{rtol} \) based on meaningful relative state error (e.g., 0.1% or 1%).
• Choose \( \mathrm{atol} \) to protect small-magnitude states from division by tiny values in the error norm.
• If a state has physical units and tight safety constraints (e.g., current, pressure, displacement), define tolerances in those physical scales.
• Inspect accepted step sizes as a diagnostic: sudden reductions often indicate rapid transients or poorly scaled states.

For LTI models with known dominant time constants, a fixed-step RK4 simulation often uses a step much smaller than the smallest dominant time constant. Adaptive methods are preferred when the dynamics are multi-scale but not stiff.

The following file implements RK2, RK4, and adaptive RK4 step-doubling for the test IVP \( y' = -2y + \sin t \), with an exact solution for error verification. It also compares against scipy.integrate.solve_ivp using the library RK45 integrator when SciPy is available.

Chapter15_Lesson2.py

# Chapter15_Lesson2.py
# Runge-Kutta Methods and Step Size Selection
# University-level demonstration: RK2, RK4, and adaptive RK4 (step-doubling)
# Optional comparison with SciPy solve_ivp (RK45) if SciPy is installed.

import math
from dataclasses import dataclass
from typing import Callable, List, Tuple

try:
    import matplotlib.pyplot as plt
except Exception:
    plt = None

try:
    from scipy.integrate import solve_ivp  # Optional library comparison
except Exception:
    solve_ivp = None


def f(t: float, y: float) -> float:
    # Test IVP: y' = -2 y + sin(t), y(0)=1
    return -2.0 * y + math.sin(t)


def y_exact(t: float) -> float:
    return (6.0 / 5.0) * math.exp(-2.0 * t) + (2.0 * math.sin(t) - math.cos(t)) / 5.0


def rk2_midpoint_step(fun: Callable[[float, float], float], t: float, y: float, h: float) -> float:
    k1 = fun(t, y)
    k2 = fun(t + 0.5 * h, y + 0.5 * h * k1)
    return y + h * k2


def rk4_step(fun: Callable[[float, float], float], t: float, y: float, h: float) -> float:
    k1 = fun(t, y)
    k2 = fun(t + 0.5 * h, y + 0.5 * h * k1)
    k3 = fun(t + 0.5 * h, y + 0.5 * h * k2)
    k4 = fun(t + h, y + h * k3)
    return y + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def integrate_fixed(
    step_fun: Callable[[Callable[[float, float], float], float, float, float], float],
    fun: Callable[[float, float], float],
    t0: float,
    tf: float,
    y0: float,
    h: float,
) -> Tuple[List[float], List[float]]:
    t = t0
    y = y0
    ts = [t]
    ys = [y]
    n_steps = int(math.ceil((tf - t0) / h))
    for _ in range(n_steps):
        h_step = min(h, tf - t)
        if h_step <= 0:
            break
        y = step_fun(fun, t, y, h_step)
        t += h_step
        ts.append(t)
        ys.append(y)
    return ts, ys


@dataclass
class AdaptiveResult:
    t: List[float]
    y: List[float]
    h_history: List[float]
    accepted: int
    rejected: int


def integrate_adaptive_rk4_stepdoubling(
    fun: Callable[[float, float], float],
    t0: float,
    tf: float,
    y0: float,
    h0: float = 0.2,
    atol: float = 1e-7,
    rtol: float = 1e-6,
    h_min: float = 1e-8,
    h_max: float = 0.5,
    safety: float = 0.9,
) -> AdaptiveResult:
    """
    Adaptive controller using RK4 with step doubling:
      y_full = RK4(h)
      y_half2 = RK4(h/2) composed twice
      error estimate (scalar) ~= |y_half2 - y_full| / (2^4 - 1) = ... / 15
    The accepted state uses Richardson-improved estimate:
      y_acc = y_half2 + (y_half2 - y_full)/15
    """
    t = t0
    y = y0
    h = h0

    ts = [t]
    ys = [y]
    hs = []
    accepted = 0
    rejected = 0

    while t < tf:
        h = min(h, tf - t)
        if h < h_min:
            raise RuntimeError("Step size fell below h_min; tolerance may be too strict.")

        y_full = rk4_step(fun, t, y, h)

        y_half = rk4_step(fun, t, y, 0.5 * h)
        y_half2 = rk4_step(fun, t + 0.5 * h, y_half, 0.5 * h)

        err_est = abs(y_half2 - y_full) / 15.0
        scale = atol + rtol * max(abs(y), abs(y_half2))
        err_norm = err_est / scale if scale > 0 else err_est

        if err_norm <= 1.0:
            # Accept
            y = y_half2 + (y_half2 - y_full) / 15.0
            t += h
            ts.append(t)
            ys.append(y)
            hs.append(h)
            accepted += 1

            # Local error estimate behaves like O(h^5) for RK4, so exponent is 1/5
            if err_norm == 0.0:
                factor = 2.0
            else:
                factor = safety * (1.0 / err_norm) ** (1.0 / 5.0)
            factor = min(2.0, max(0.2, factor))
            h = min(h_max, factor * h)
        else:
            # Reject and reduce step
            rejected += 1
            factor = safety * (1.0 / max(err_norm, 1e-16)) ** (1.0 / 5.0)
            factor = min(1.0, max(0.1, factor))
            h = max(h_min, factor * h)

    return AdaptiveResult(ts, ys, hs, accepted, rejected)


def max_error(ts: List[float], ys: List[float]) -> float:
    return max(abs(y - y_exact(t)) for t, y in zip(ts, ys))


def main() -> None:
    t0, tf, y0 = 0.0, 10.0, 1.0

    # Fixed-step runs
    h = 0.1
    t_rk2, y_rk2 = integrate_fixed(rk2_midpoint_step, f, t0, tf, y0, h)
    t_rk4, y_rk4 = integrate_fixed(rk4_step, f, t0, tf, y0, h)

    # Adaptive RK4 (step-doubling)
    adaptive = integrate_adaptive_rk4_stepdoubling(f, t0, tf, y0, h0=0.2, atol=1e-8, rtol=1e-6)

    print("=== Fixed-step results ===")
    print(f"RK2 midpoint, h={h:.3f}, steps={len(t_rk2)-1}, max error={max_error(t_rk2, y_rk2):.6e}")
    print(f"RK4,         h={h:.3f}, steps={len(t_rk4)-1}, max error={max_error(t_rk4, y_rk4):.6e}")

    print("\n=== Adaptive RK4 (step-doubling) ===")
    print(f"Accepted steps: {adaptive.accepted}")
    print(f"Rejected steps: {adaptive.rejected}")
    print(f"Max error:      {max_error(adaptive.t, adaptive.y):.6e}")
    if adaptive.h_history:
        print(f"h range:        [{min(adaptive.h_history):.3e}, {max(adaptive.h_history):.3e}]")

    if solve_ivp is not None:
        sol = solve_ivp(
            lambda t, y: [-2.0 * y[0] + math.sin(t)],
            (t0, tf),
            [y0],
            method="RK45",
            atol=1e-8,
            rtol=1e-6,
            dense_output=True,
        )
        err_scipy = 0.0
        for t in adaptive.t:
            yv = float(sol.sol(t)[0])
            err_scipy = max(err_scipy, abs(yv - y_exact(t)))
        print("\n=== SciPy solve_ivp (RK45) comparison ===")
        print(f"Status: {sol.message}")
        print(f"Internal steps: {len(sol.t)-1}")
        print(f"Max error on adaptive grid: {err_scipy:.6e}")
    else:
        print("\nSciPy not installed; skipping library comparison.")

    if plt is not None:
        t_dense = [i * (tf / 1000.0) for i in range(1001)]
        y_dense = [y_exact(t) for t in t_dense]

        plt.figure()
        plt.plot(t_dense, y_dense, label="Exact")
        plt.plot(t_rk2, y_rk2, ".", label="RK2 (h=0.1)")
        plt.plot(t_rk4, y_rk4, "-", label="RK4 (h=0.1)")
        plt.plot(adaptive.t, adaptive.y, "o", markersize=3, label="Adaptive RK4")
        plt.xlabel("t")
        plt.ylabel("y(t)")
        plt.title("Chapter15_Lesson2: RK Methods")
        plt.legend()
        plt.grid(True)

        plt.figure()
        if adaptive.h_history:
            plt.plot(adaptive.t[1:], adaptive.h_history)
            plt.xlabel("t")
            plt.ylabel("Accepted step size h")
            plt.title("Adaptive step-size history")
            plt.grid(True)
        plt.show()


if __name__ == "__main__":
    main()

      

Library note: In Python, the most common high-level API is scipy.integrate.solve_ivp, with methods such as RK45 and DOP853 for non-stiff problems.

This C++17 code implements fixed-step RK4 and adaptive RK4 step-doubling. The file is self-contained and compiles without external dependencies.

Chapter15_Lesson2.cpp

// Chapter15_Lesson2.cpp
// Runge-Kutta Methods and Step Size Selection
// C++17 example: RK4 + adaptive step-doubling controller.
// Related library (not required to compile this file): Boost.Odeint.

#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <utility>
#include <vector>

struct AdaptiveResult {
    std::vector<double> t;
    std::vector<double> y;
    std::vector<double> h_history;
    int accepted = 0;
    int rejected = 0;
};

double f(double t, double y) {
    return -2.0 * y + std::sin(t);
}

double y_exact(double t) {
    return (6.0 / 5.0) * std::exp(-2.0 * t) + (2.0 * std::sin(t) - std::cos(t)) / 5.0;
}

double rk4_step(double (*fun)(double, double), double t, double y, double h) {
    const double k1 = fun(t, y);
    const double k2 = fun(t + 0.5 * h, y + 0.5 * h * k1);
    const double k3 = fun(t + 0.5 * h, y + 0.5 * h * k2);
    const double k4 = fun(t + h, y + h * k3);
    return y + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

std::pair<std::vector<double>, std::vector<double>> integrate_fixed_rk4(
    double t0, double tf, double y0, double h
) {
    std::vector<double> ts{t0};
    std::vector<double> ys{y0};
    double t = t0;
    double y = y0;

    while ((tf - t) > 1e-12) {
        const double h_step = std::min(h, tf - t);
        y = rk4_step(f, t, y, h_step);
        t += h_step;
        ts.push_back(t);
        ys.push_back(y);
    }
    return {ts, ys};
}

AdaptiveResult integrate_adaptive_rk4_stepdoubling(
    double t0, double tf, double y0,
    double h0 = 0.2, double atol = 1e-8, double rtol = 1e-6,
    double h_min = 1e-8, double h_max = 0.5, double safety = 0.9
) {
    AdaptiveResult out;
    out.t.push_back(t0);
    out.y.push_back(y0);

    double t = t0;
    double y = y0;
    double h = h0;

    while ((tf - t) > 1e-12) {
        h = std::min(h, tf - t);
        if (h < h_min) {
            throw std::runtime_error("Step size below h_min");
        }

        const double y_full = rk4_step(f, t, y, h);
        const double y_half = rk4_step(f, t, y, 0.5 * h);
        const double y_half2 = rk4_step(f, t + 0.5 * h, y_half, 0.5 * h);

        const double err_est = std::abs(y_half2 - y_full) / 15.0;
        const double scale = atol + rtol * std::max(std::abs(y), std::abs(y_half2));
        const double err_norm = (scale > 0.0) ? (err_est / scale) : err_est;

        if (err_norm <= 1.0) {
            y = y_half2 + (y_half2 - y_full) / 15.0;
            t += h;
            out.t.push_back(t);
            out.y.push_back(y);
            out.h_history.push_back(h);
            out.accepted++;

            double factor;
            if (err_norm == 0.0) {
                factor = 2.0;
            } else {
                factor = safety * std::pow(1.0 / err_norm, 1.0 / 5.0);
            }
            factor = std::min(2.0, std::max(0.2, factor));
            h = std::min(h_max, factor * h);
        } else {
            out.rejected++;
            double factor = safety * std::pow(1.0 / std::max(err_norm, 1e-16), 1.0 / 5.0);
            factor = std::min(1.0, std::max(0.1, factor));
            h = std::max(h_min, factor * h);
        }
    }

    return out;
}

double max_error(const std::vector<double>& ts, const std::vector<double>& ys) {
    double emax = 0.0;
    for (std::size_t i = 0; i < ts.size(); ++i) {
        emax = std::max(emax, std::abs(ys[i] - y_exact(ts[i])));
    }
    return emax;
}

int main() {
    const double t0 = 0.0, tf = 10.0, y0 = 1.0;

    auto [t_fixed, y_fixed] = integrate_fixed_rk4(t0, tf, y0, 0.1);
    AdaptiveResult adaptive = integrate_adaptive_rk4_stepdoubling(t0, tf, y0);

    std::cout << std::scientific << std::setprecision(6);
    std::cout << "Chapter15_Lesson2.cpp\n";
    std::cout << "Fixed RK4 (h=0.1): steps = " << (t_fixed.size() - 1)
              << ", max error = " << max_error(t_fixed, y_fixed) << "\n";

    std::cout << "Adaptive RK4 step-doubling: accepted = " << adaptive.accepted
              << ", rejected = " << adaptive.rejected
              << ", max error = " << max_error(adaptive.t, adaptive.y) << "\n";

    if (!adaptive.h_history.empty()) {
        double hmin = adaptive.h_history.front();
        double hmax = adaptive.h_history.front();
        for (double h : adaptive.h_history) {
            hmin = std::min(hmin, h);
            hmax = std::max(hmax, h);
        }
        std::cout << "h range = [" << hmin << ", " << hmax << "]\n";
    }

    std::cout << "\nLibrary note: Boost.Odeint offers explicit_rk4, runge_kutta_cash_karp54, "
                 "and controlled steppers for adaptive integration.\n";
    return 0;
}

      

Library note: In production C++ workflows, Boost.Odeint provides explicit RK steppers and controlled adaptive steppers (e.g., Cash–Karp 5(4)).

The Java implementation mirrors the same algorithmic structure, which is useful for students building simulation tools in JVM-based environments.

Chapter15_Lesson2.java

// Chapter15_Lesson2.java
// Runge-Kutta Methods and Step Size Selection
// Java example: RK4 + adaptive step-doubling controller.
// Related library: Apache Commons Math (DormandPrince54Integrator).

import java.util.ArrayList;
import java.util.List;

public class Chapter15_Lesson2 {

    static class AdaptiveResult {
        List<Double> t = new ArrayList<>();
        List<Double> y = new ArrayList<>();
        List<Double> hHistory = new ArrayList<>();
        int accepted = 0;
        int rejected = 0;
    }

    static double f(double t, double y) {
        return -2.0 * y + Math.sin(t);
    }

    static double yExact(double t) {
        return (6.0 / 5.0) * Math.exp(-2.0 * t) + (2.0 * Math.sin(t) - Math.cos(t)) / 5.0;
    }

    static double rk4Step(double t, double y, double h) {
        double k1 = f(t, y);
        double k2 = f(t + 0.5 * h, y + 0.5 * h * k1);
        double k3 = f(t + 0.5 * h, y + 0.5 * h * k2);
        double k4 = f(t + h, y + h * k3);
        return y + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
    }

    static double[] integrateFixedRK4(double t0, double tf, double y0, double h) {
        int n = (int)Math.ceil((tf - t0) / h);
        double[] out = new double[n + 1];  // store y values only (grid is uniform)
        out[0] = y0;
        double t = t0;
        double y = y0;
        for (int i = 1; i <= n; i++) {
            double hs = Math.min(h, tf - t);
            y = rk4Step(t, y, hs);
            t += hs;
            out[i] = y;
        }
        return out;
    }

    static AdaptiveResult integrateAdaptiveRK4StepDoubling(
            double t0, double tf, double y0,
            double h0, double atol, double rtol,
            double hMin, double hMax, double safety) {

        AdaptiveResult res = new AdaptiveResult();
        res.t.add(t0);
        res.y.add(y0);

        double t = t0;
        double y = y0;
        double h = h0;

        while (t < tf) {
            h = Math.min(h, tf - t);
            if (h < hMin) {
                throw new RuntimeException("Step size below hMin");
            }

            double yFull = rk4Step(t, y, h);
            double yHalf = rk4Step(t, y, 0.5 * h);
            double yHalf2 = rk4Step(t + 0.5 * h, yHalf, 0.5 * h);

            double errEst = Math.abs(yHalf2 - yFull) / 15.0;
            double scale = atol + rtol * Math.max(Math.abs(y), Math.abs(yHalf2));
            double errNorm = (scale > 0.0) ? errEst / scale : errEst;

            if (errNorm <= 1.0) {
                y = yHalf2 + (yHalf2 - yFull) / 15.0;
                t += h;
                res.t.add(t);
                res.y.add(y);
                res.hHistory.add(h);
                res.accepted++;

                double factor;
                if (errNorm == 0.0) {
                    factor = 2.0;
                } else {
                    factor = safety * Math.pow(1.0 / errNorm, 1.0 / 5.0);
                }
                factor = Math.min(2.0, Math.max(0.2, factor));
                h = Math.min(hMax, factor * h);
            } else {
                res.rejected++;
                double factor = safety * Math.pow(1.0 / Math.max(errNorm, 1e-16), 1.0 / 5.0);
                factor = Math.min(1.0, Math.max(0.1, factor));
                h = Math.max(hMin, factor * h);
            }
        }

        return res;
    }

    static double maxErrorUniformGrid(double t0, double h, double[] ys) {
        double emax = 0.0;
        for (int i = 0; i < ys.length; i++) {
            double t = t0 + i * h;
            emax = Math.max(emax, Math.abs(ys[i] - yExact(t)));
        }
        return emax;
    }

    static double maxErrorAdaptive(AdaptiveResult res) {
        double emax = 0.0;
        for (int i = 0; i < res.t.size(); i++) {
            emax = Math.max(emax, Math.abs(res.y.get(i) - yExact(res.t.get(i))));
        }
        return emax;
    }

    public static void main(String[] args) {
        double t0 = 0.0, tf = 10.0, y0 = 1.0, h = 0.1;

        double[] yFixed = integrateFixedRK4(t0, tf, y0, h);
        AdaptiveResult adaptive = integrateAdaptiveRK4StepDoubling(
                t0, tf, y0, 0.2, 1e-8, 1e-6, 1e-8, 0.5, 0.9);

        System.out.println("Chapter15_Lesson2.java");
        System.out.printf("Fixed RK4 (h=%.3f): steps=%d, max error=%.6e%n",
                h, yFixed.length - 1, maxErrorUniformGrid(t0, h, yFixed));

        System.out.printf("Adaptive RK4 step-doubling: accepted=%d, rejected=%d, max error=%.6e%n",
                adaptive.accepted, adaptive.rejected, maxErrorAdaptive(adaptive));

        if (!adaptive.hHistory.isEmpty()) {
            double hMin = adaptive.hHistory.get(0);
            double hMax = adaptive.hHistory.get(0);
            for (double hs : adaptive.hHistory) {
                hMin = Math.min(hMin, hs);
                hMax = Math.max(hMax, hs);
            }
            System.out.printf("h range=[%.3e, %.3e]%n", hMin, hMax);
        }

        System.out.println("Library note: Apache Commons Math provides "
                + "DormandPrince54Integrator and ClassicalRungeKuttaIntegrator.");
    }
}

      

Library note: Apache Commons Math includes ClassicalRungeKuttaIntegrator and DormandPrince54Integrator for fixed and adaptive explicit integration.

The MATLAB file below implements RK2, RK4, and adaptive RK4 step-doubling, and compares the result with ode45 (Dormand–Prince family). A short Simulink setup note is included in the comments.

Chapter15_Lesson2.m

% Chapter15_Lesson2.m
% Runge-Kutta Methods and Step Size Selection
% MATLAB script: RK2, RK4, adaptive RK4 (step-doubling), and ode45 comparison.

clear; clc;

f = @(t,y) -2*y + sin(t);
yExact = @(t) (6/5)*exp(-2*t) + (2*sin(t) - cos(t))/5;

t0 = 0; tf = 10; y0 = 1;

% Fixed-step RK2 and RK4
h = 0.1;
[t2, y2] = integrateFixed(@rk2midpointStep, f, t0, tf, y0, h);
[t4, y4] = integrateFixed(@rk4Step,       f, t0, tf, y0, h);

% Adaptive RK4 with step-doubling
res = integrateAdaptiveRK4StepDoubling(f, t0, tf, y0, 0.2, 1e-8, 1e-6, 1e-8, 0.5, 0.9);

% MATLAB built-in adaptive RK pair (Dormand-Prince)
opts = odeset('RelTol',1e-6,'AbsTol',1e-8);
[t45, y45] = ode45(f, [t0 tf], y0, opts);

fprintf('Chapter15_Lesson2.m\n');
fprintf('RK2 midpoint (h=%.3f), max error = %.6e\n', h, max(abs(y2 - yExact(t2))));
fprintf('RK4         (h=%.3f), max error = %.6e\n', h, max(abs(y4 - yExact(t4))));
fprintf('Adaptive RK4 step-doubling: accepted=%d, rejected=%d, max error=%.6e\n', ...
    res.accepted, res.rejected, max(abs(res.y - yExact(res.t))));
fprintf('ode45 internal steps = %d, max error = %.6e\n', length(t45)-1, max(abs(y45 - yExact(t45))));

figure;
tt = linspace(t0, tf, 1000);
plot(tt, yExact(tt), 'LineWidth', 1.5); hold on;
plot(t2, y2, '.');
plot(t4, y4, '-');
plot(res.t, res.y, 'o', 'MarkerSize', 3);
legend('Exact','RK2','RK4','Adaptive RK4');
xlabel('t'); ylabel('y(t)'); grid on;
title('Chapter15\_Lesson2: RK Methods');

figure;
plot(res.t(2:end), res.hHistory, '-o', 'MarkerSize', 3);
xlabel('t'); ylabel('Accepted step size h'); grid on;
title('Adaptive step-size history');

% ---- Simulink note ----
% In Simulink, build the same scalar ODE using:
%   Integrator block + Sum block + Sine Wave block + Gain(-2)
% and compare fixed-step "ode4 (Runge-Kutta)" vs variable-step "ode45".

%% Local functions
function yNext = rk2midpointStep(fun, t, y, h)
    k1 = fun(t, y);
    k2 = fun(t + 0.5*h, y + 0.5*h*k1);
    yNext = y + h*k2;
end

function yNext = rk4Step(fun, t, y, h)
    k1 = fun(t, y);
    k2 = fun(t + 0.5*h, y + 0.5*h*k1);
    k3 = fun(t + 0.5*h, y + 0.5*h*k2);
    k4 = fun(t + h, y + h*k3);
    yNext = y + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
end

function [ts, ys] = integrateFixed(stepFun, fun, t0, tf, y0, h)
    ts = t0:h:tf;
    if ts(end) < tf
        ts = [ts tf];
    end
    ys = zeros(size(ts));
    ys(1) = y0;
    for i = 1:length(ts)-1
        hs = ts(i+1) - ts(i);
        ys(i+1) = stepFun(fun, ts(i), ys(i), hs);
    end
end

function res = integrateAdaptiveRK4StepDoubling(fun, t0, tf, y0, h0, atol, rtol, hMin, hMax, safety)
    t = t0; y = y0; h = h0;
    ts = t; ys = y; hHistory = [];
    accepted = 0; rejected = 0;

    while t < tf
        h = min(h, tf - t);
        if h < hMin
            error('Step size below hMin');
        end

        yFull = rk4Step(fun, t, y, h);
        yHalf = rk4Step(fun, t, y, 0.5*h);
        yHalf2 = rk4Step(fun, t + 0.5*h, yHalf, 0.5*h);

        errEst = abs(yHalf2 - yFull) / 15;
        scale = atol + rtol * max(abs(y), abs(yHalf2));
        errNorm = errEst / scale;

        if errNorm <= 1
            y = yHalf2 + (yHalf2 - yFull)/15;
            t = t + h;
            ts(end+1,1) = t; %#ok<AGROW>
            ys(end+1,1) = y; %#ok<AGROW>
            hHistory(end+1,1) = h; %#ok<AGROW>
            accepted = accepted + 1;

            if errNorm == 0
                factor = 2;
            else
                factor = safety * (1/errNorm)^(1/5);
            end
            factor = min(2, max(0.2, factor));
            h = min(hMax, factor*h);
        else
            rejected = rejected + 1;
            factor = safety * (1/max(errNorm,1e-16))^(1/5);
            factor = min(1, max(0.1, factor));
            h = max(hMin, factor*h);
        end
    end

    res.t = ts;
    res.y = ys;
    res.hHistory = hHistory;
    res.accepted = accepted;
    res.rejected = rejected;
end

      

Simulink note: Use the model-level solver settings to compare fixed-step ode4 (Runge-Kutta) and variable-step ode45. This is a direct way to visualize how solver choice changes step density over transients.

The following Wolfram Language code (packaged in the course as Chapter15_Lesson2.nb) implements the same RK methods and compares them with NDSolve using an explicit Runge–Kutta method.

Chapter15_Lesson2.nb

(* Chapter15_Lesson2.nb *)
(* Runge-Kutta Methods and Step Size Selection *)
(* Wolfram Language code (saved with .nb filename for course packaging) *)

ClearAll["Global`*"];

f[t_, y_] := -2 y + Sin[t];
yExact[t_] := (6/5) Exp[-2 t] + (2 Sin[t] - Cos[t])/5;

rk2MidpointStep[fun_, t_, y_, h_] := Module[{k1, k2},
  k1 = fun[t, y];
  k2 = fun[t + h/2, y + (h/2) k1];
  y + h k2
];

rk4Step[fun_, t_, y_, h_] := Module[{k1, k2, k3, k4},
  k1 = fun[t, y];
  k2 = fun[t + h/2, y + (h/2) k1];
  k3 = fun[t + h/2, y + (h/2) k2];
  k4 = fun[t + h, y + h k3];
  y + (h/6) (k1 + 2 k2 + 2 k3 + k4)
];

integrateFixed[stepFun_, fun_, {t0_, tf_}, y0_, h_] := Module[
  {t = t0, y = y0, ts = {t0}, ys = {y0}, hs},
  While[t < tf,
    hs = Min[h, tf - t];
    y = stepFun[fun, t, y, hs];
    t = t + hs;
    AppendTo[ts, t];
    AppendTo[ys, y];
  ];
  <|"t" -> ts, "y" -> ys|>
];

integrateAdaptiveRK4StepDoubling[fun_, {t0_, tf_}, y0_, h0_: 0.2, atol_: 1.*^-8, rtol_: 1.*^-6] := Module[
  {t = t0, y = y0, h = h0, hMin = 1.*^-8, hMax = 0.5, safety = 0.9,
   ts = {t0}, ys = {y0}, hHistory = {}, accepted = 0, rejected = 0,
   yFull, yHalf, yHalf2, errEst, scale, errNorm, factor},
  While[t < tf,
    h = Min[h, tf - t];
    If[h < hMin, Throw["Step size below hMin"]];
    yFull = rk4Step[fun, t, y, h];
    yHalf = rk4Step[fun, t, y, h/2];
    yHalf2 = rk4Step[fun, t + h/2, yHalf, h/2];

    errEst = Abs[yHalf2 - yFull]/15;
    scale = atol + rtol Max[Abs[y], Abs[yHalf2]];
    errNorm = errEst/scale;

    If[errNorm <= 1,
      y = yHalf2 + (yHalf2 - yFull)/15;
      t = t + h;
      AppendTo[ts, t];
      AppendTo[ys, y];
      AppendTo[hHistory, h];
      accepted++;
      factor = If[errNorm == 0, 2, safety (1/errNorm)^(1/5)];
      factor = Min[2, Max[0.2, factor]];
      h = Min[hMax, factor h];,
      rejected++;
      factor = safety (1/Max[errNorm, 10^-16])^(1/5);
      factor = Min[1, Max[0.1, factor]];
      h = Max[hMin, factor h];
    ];
  ];
  <|"t" -> ts, "y" -> ys, "hHistory" -> hHistory, "accepted" -> accepted, "rejected" -> rejected|>
];

t0 = 0; tf = 10; y0 = 1; h = 0.1;
solRK2 = integrateFixed[rk2MidpointStep, f, {t0, tf}, y0, h];
solRK4 = integrateFixed[rk4Step, f, {t0, tf}, y0, h];
solAdaptive = integrateAdaptiveRK4StepDoubling[f, {t0, tf}, y0];

maxErr[data_] := Max[Abs[data["y"] - (yExact /@ data["t"])]];
Print["Chapter15_Lesson2.nb"];
Print["RK2 midpoint max error = ", ScientificForm[maxErr[solRK2], 6]];
Print["RK4 max error = ", ScientificForm[maxErr[solRK4], 6]];
Print["Adaptive RK4 max error = ", ScientificForm[maxErr[solAdaptive], 6]];
Print["Accepted steps = ", solAdaptive["accepted"], ", Rejected steps = ", solAdaptive["rejected"]];

(* Library-style comparison in Mathematica *)
nd = NDSolveValue[{y'[t] == f[t, y[t]], y[0] == 1}, y, {t, 0, 10},
  Method -> {"ExplicitRungeKutta"}];

ListLinePlot[
 {
  Transpose[{solRK4["t"], solRK4["y"]}],
  Table[{tt, yExact[tt]}, {tt, 0, 10, 0.01}],
  Table[{tt, nd[tt]}, {tt, 0, 10, 0.01}]
 },
 PlotLegends -> {"RK4", "Exact", "NDSolve ExplicitRungeKutta"},
 GridLines -> Automatic
]

      

Library note: Mathematica's NDSolve is the standard high-level solver interface; specifying Method -> {"ExplicitRungeKutta"} is useful for controlled comparisons in teaching.

Let the exact one-step result be \( y^* \). For RK4, the local one-step truncation error has the form

\[ y^* - y_{\text{full}} = C h^5 + O(h^6) \]

For two half-steps, each contributes an error proportional to \( (h/2)^5 \), so the total is

\[ y^* - y_{\text{half2}} = 2 C \left(\frac{h}{2}\right)^5 + O(h^6) = \frac{C h^5}{16} + O(h^6) \]

Subtracting gives

\[ y_{\text{half2}} - y_{\text{full}} = \frac{15}{16} C h^5 + O(h^6) \]

Therefore

\[ C h^5 = \frac{16}{15}\left(y_{\text{half2}} - y_{\text{full}}\right) + O(h^6) \]

and the more accurate (Richardson-improved) accepted value is

\[ y_{\text{acc}} = y_{\text{half2}} + \frac{y_{\text{half2}} - y_{\text{full}}}{15} \]

which cancels the leading \( O(h^5) \) term.

Problem 1 (Midpoint RK2 order proof): Starting from \( y' = f(t,y) \), prove that the midpoint method has local truncation error \( O(h^3) \).

Solution: Use the expansion \( k_2 = f_n + \frac{h}{2}(f_t + f_y f)_n + O(h^2) \). Then

\[ y_{n+1} = y_n + h k_2 = y_n + h f_n + \frac{h^2}{2}(f_t + f_y f)_n + O(h^3) \]

Since \( y'' = f_t + f_y f \) on the exact trajectory, this matches the exact Taylor series through the \( h^2 \) term. Hence the local truncation error is \( O(h^3) \), and the global error is \( O(h^2) \).

Problem 2 (Observed order test): Let \( Y_h(T) \), \( Y_{h/2}(T) \), and \( Y_{h/4}(T) \) be numerical approximations at final time \( T \) from a method with global error \( K h^p \). Show that the observed order estimate

\[ p_{\mathrm{obs}} = \log_2\!\left( \frac{|Y_h(T)-Y_{h/2}(T)|}{|Y_{h/2}(T)-Y_{h/4}(T)|} \right) \]

converges to \( p \).

Solution: Under the asymptotic model,

\[ Y_h - Y_{h/2} = K h^p\left(1-2^{-p}\right), \qquad Y_{h/2} - Y_{h/4} = K (h/2)^p\left(1-2^{-p}\right) \]

so their ratio tends to \( 2^p \). Taking \( \log_2 \) yields \( p_{\mathrm{obs}} \to p \).

Problem 3 (Adaptive controller exponent): Suppose the estimated local error behaves as \( e(h)=C h^q \). Derive the formula for the next step size that targets \( e(h_{\text{new}})=1 \) (after normalization).

Solution: Set \( C h_{\text{new}}^q = 1 \) and use \( C = e(h)/h^q \). Then

\[ h_{\text{new}} = h\left(\frac{1}{e(h)}\right)^{1/q} \]

In practice, include a safety factor: \( h_{\text{new}} = s h (1/e)^{1/q} \). For RK4 step-doubling, \( q=5 \).

Problem 4 (Scalar test equation and RK4 polynomial): For the test equation \( y' = \lambda y \), show that one RK4 step gives \( y_{n+1} = R(z) y_n \) with \( z = h\lambda \) and

\[ R(z) = 1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24} \]

Solution: Substitute \( f(t,y)=\lambda y \) into the RK4 stages. Each stage becomes a polynomial multiple of \( y_n \), and after summation the update factor is exactly the fourth Taylor polynomial of \( e^z \). This confirms fourth-order matching for the linear test equation and also provides intuition for numerical amplification.

Problem 5 (Tolerance scaling for multi-state models): A system has two states: displacement \( x \) in meters and current \( i \) in amperes. Explain why using one absolute tolerance for both states can be problematic, and propose a better strategy.

Solution: The state magnitudes and physical units differ, so a single absolute tolerance may over-penalize one state and under-penalize the other. Use a componentwise scale \( \mathrm{sc}_i = \mathrm{atol}_i + \mathrm{rtol}\max(|y_i|,|\hat{y}_i|) \) with \( \mathrm{atol}_i \) chosen per physical variable (e.g., micrometers for displacement and milliamps for current). This makes the adaptive controller physically meaningful.

We derived RK2 and RK4 from Taylor-series ideas, introduced order conditions, and proved the local-to-global error relationship using a discrete Gronwall argument. We then developed a practical adaptive step-size mechanism based on RK4 step-doubling, including the 1/15 Richardson correction and a safety-factor controller. The provided Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica implementations give a consistent computational template for non-stiff IVPs in system dynamics.