In system dynamics, most continuous-time models are written as an initial value problem in first-order form:

\[ \dot{\mathbf{x}}(t) = \mathbf{f}(t,\mathbf{x}(t),\mathbf{u}(t)), \qquad \mathbf{x}(t_0)=\mathbf{x}_0 \]

Here \( \mathbf{x}(t)\in\mathbb{R}^n \) is the state vector, \( \mathbf{u}(t) \) is the known input, and \( \mathbf{f} \) encodes the physics of the system. For example, a mass–spring–damper system \( m\ddot{q}+c\dot{q}+kq=F(t) \) becomes first-order by defining \( x_1=q \) and \( x_2=\dot{q} \):

\[ \dot{\mathbf{x}} = \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \begin{bmatrix} x_2 \\ -\dfrac{k}{m}x_1 - \dfrac{c}{m}x_2 + \dfrac{1}{m}F(t) \end{bmatrix}, \qquad \mathbf{x}(t_0)= \begin{bmatrix} q(t_0) \\ \dot{q}(t_0) \end{bmatrix}. \]

Numerical integration is required whenever a closed-form solution is unavailable, too expensive, or not useful for simulation workflows (parameter sweeps, control testing, digital twins, etc.).

A standard well-posedness condition is local Lipschitz continuity in \( \mathbf{x} \). If \( \mathbf{f} \) is continuous in \( t \) and locally Lipschitz in \( \mathbf{x} \), then a unique local solution exists for the IVP. This theoretical result is important because numerical methods approximate a solution only when the underlying IVP itself is well-defined.

flowchart TD
  A["Physical laws (Newton, KVL, energy balance)"] --> B["Build ODE model"]
  B --> C["Convert to first-order state form"]
  C --> D["Set initial state x(t0) and input u(t)"]
  D --> E["Choose step size h and method"]
  E --> F["Propagate state on grid t0, t1, ..., tN"]
  F --> G["Check error and stability"]
  G --> H["Accept results or refine h"]
        

We discretize time on a uniform grid:

\[ t_k = t_0 + kh, \qquad k=0,1,\dots,N, \qquad h=\dfrac{T-t_0}{N} \]

Let \( \mathbf{x}_k \approx \mathbf{x}(t_k) \) denote the numerical approximation. A one-step method has the general form

\[ \mathbf{x}_{k+1} = \mathbf{x}_k + h\,\boldsymbol{\Phi}(t_k,\mathbf{x}_k,h) \]

where \( \boldsymbol{\Phi} \) is the method-specific slope model. A method is consistent if, informally, its one-step formula reproduces the differential equation as \( h \to 0 \). For smooth solutions, consistency is established by Taylor expansion.

We will use two error notions throughout:

• Local truncation error (LTE): the defect of one numerical step when the exact state at \( t_k \) is used as input.
• Global error: the accumulated error at grid nodes after many steps.

In this lesson, we will prove that explicit Euler is first-order accurate globally \( O(h) \), while improved Euler (Heun) is second-order globally \( O(h^2) \).

Starting from Taylor expansion of the exact solution at \( t_k \):

\[ \mathbf{x}(t_{k+1}) = \mathbf{x}(t_k) + h\dot{\mathbf{x}}(t_k) + \dfrac{h^2}{2}\ddot{\mathbf{x}}(\xi_k), \qquad \xi_k \in (t_k,t_{k+1}) \]

and using \( \dot{\mathbf{x}}(t_k)=\mathbf{f}(t_k,\mathbf{x}(t_k),\mathbf{u}(t_k)) \), the first-order truncation of the expansion gives the explicit Euler update:

\[ \mathbf{x}_{k+1} = \mathbf{x}_k + h\,\mathbf{f}(t_k,\mathbf{x}_k,\mathbf{u}(t_k)) \]

\( \text{Local truncation error} \) is obtained by inserting the exact solution into the scheme:

\[ \boldsymbol{\eta}_{k+1}^{(E)} = \mathbf{x}(t_{k+1}) - \mathbf{x}(t_k) - h\,\mathbf{f}(t_k,\mathbf{x}(t_k),\mathbf{u}(t_k)) = \dfrac{h^2}{2}\ddot{\mathbf{x}}(\xi_k) \]

Therefore, if \( \ddot{\mathbf{x}} \) is bounded on the interval, then \( \|\boldsymbol{\eta}_{k+1}^{(E)}\| \le C h^2 \), so Euler has local order 2 and global order 1.

To show the global error bound, define \( \mathbf{e}_k = \mathbf{x}(t_k)-\mathbf{x}_k \). Assume \( \mathbf{f} \) is Lipschitz in \( \mathbf{x} \) with constant \( L \). Subtracting the Euler recursion from the exact one-step identity gives

\[ \mathbf{e}_{k+1} = \mathbf{e}_k + h\!\left[\mathbf{f}(t_k,\mathbf{x}(t_k),\mathbf{u}(t_k))-\mathbf{f}(t_k,\mathbf{x}_k,\mathbf{u}(t_k))\right] + \boldsymbol{\eta}_{k+1}^{(E)} \]

Taking norms and applying the Lipschitz property:

\[ \|\mathbf{e}_{k+1}\| \le (1+hL)\|\mathbf{e}_k\| + Ch^2 \]

Since \( \mathbf{e}_0=\mathbf{0} \), repeated substitution and the discrete Grönwall inequality yield

\[ \|\mathbf{e}_k\| \le \dfrac{C}{L}\left(e^{L(t_k-t_0)}-1\right)h \]

Hence Euler is globally first-order accurate: \( \max_k \|\mathbf{e}_k\| = O(h) \).

The improved Euler method (also called Heun’s method or explicit trapezoidal method) averages two slope estimates:

1. Predictor: forward Euler estimate of the next state.
2. Corrector: trapezoidal average of the slope at the start and the predicted end.

\[ \mathbf{x}_{k+1}^{(p)} = \mathbf{x}_k + h\,\mathbf{f}(t_k,\mathbf{x}_k,\mathbf{u}(t_k)) \]

\[ \mathbf{x}_{k+1} = \mathbf{x}_k + \dfrac{h}{2}\Big( \mathbf{f}(t_k,\mathbf{x}_k,\mathbf{u}(t_k)) + \mathbf{f}(t_{k+1},\mathbf{x}_{k+1}^{(p)},\mathbf{u}(t_{k+1})) \Big) \]

For a smooth scalar solution \( x(t) \), expand the exact solution and the right-hand side to compare terms. The exact Taylor series is

\[ x(t_{k+1}) = x(t_k) + h x'(t_k) + \dfrac{h^2}{2}x''(t_k) + O(h^3) \]

The predictor gives \( x_{k+1}^{(p)} = x(t_k) + h x'(t_k) + O(h^2) \) when exact input is used. Then

\[ f(t_{k+1},x_{k+1}^{(p)}) = x'(t_k) + h x''(t_k) + O(h^2) \]

Substituting into the corrector:

\[ x(t_k) + \dfrac{h}{2}\left(x'(t_k) + x'(t_k)+h x''(t_k)\right) + O(h^3) = x(t_k)+h x'(t_k)+\dfrac{h^2}{2}x''(t_k)+O(h^3) \]

which matches the exact expansion up to \( h^2 \). Thus the improved Euler method has \( \text{LTE}=O(h^3) \) and \( \text{global error}=O(h^2) \).

In system dynamics, this extra order often produces significantly better simulation accuracy at the same step size, especially for smooth trajectories.

Accuracy alone is not enough: a method can be accurate in theory but unstable if the step size is too large. The standard test problem for absolute stability is

\[ y' = \lambda y, \qquad \operatorname{Re}(\lambda) < 0 \]

Let \( z=h\lambda \). For Euler:

\[ y_{k+1} = (1+z)y_k, \qquad R_E(z)=1+z \]

Absolute stability requires the amplification factor to decay:

\[ |R_E(z)| = |1+z| < 1 \]

For negative real \( \lambda \) (pure decay mode), this becomes

\[ 0 < -h\lambda < 2 \]

For improved Euler (Heun), applying the method to \( y'=\lambda y \) yields

\[ y_{k+1} = \left(1+z+\dfrac{z^2}{2}\right)y_k, \qquad R_H(z)=1+z+\dfrac{z^2}{2} \]

and absolute stability is

\[ \left|1+z+\dfrac{z^2}{2}\right| < 1 \]

Along the negative real axis, the same interval endpoint appears:

\[ 0 < -h\lambda < 2 \]

but the two-dimensional stability region in the complex plane is different. This matters for oscillatory modes (complex eigenvalues), which frequently arise in mechanical and electrical systems.

For linear systems \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \), a practical rule is to choose \( h \) small enough relative to the fastest mode (largest negative real part or highest oscillation frequency) so that all critical eigenvalues mapped by \( z_i=h\lambda_i(\mathbf{A}) \) lie inside the method’s stability region.

flowchart TD
  S["Given IVP and model dynamics"] --> H["Choose trial step size h"]
  H --> A["Run Euler or Improved Euler"]
  A --> B["Check growth/oscillation of numerical solution"]
  B --> C{"Stable and accurate?"}
  C -->|No: unstable| D["Reduce h \n(stability issue)"]
  C -->|No: stable but inaccurate| E["Reduce h or use \nhigher-order method"]
  C -->|Yes| F["Accept simulation output"]
  D --> A
  E --> A
        

Scalar IVP. For \( y'=f(t,y) \), with \( y(t_0)=y_0 \):

\[ \text{Euler:}\qquad y_{k+1}=y_k+h f(t_k,y_k) \]

\[ \text{Improved Euler:}\qquad \begin{aligned} y_{k+1}^{(p)} &= y_k+h f(t_k,y_k) \\ y_{k+1} &= y_k+\dfrac{h}{2}\left(f(t_k,y_k)+f(t_{k+1},y_{k+1}^{(p)})\right) \end{aligned} \]

Vector IVP. For \( \dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x},\mathbf{u}) \):

\[ \text{Euler:}\qquad \mathbf{x}_{k+1}=\mathbf{x}_k+h\,\mathbf{f}(t_k,\mathbf{x}_k,\mathbf{u}_k) \]

\[ \text{Improved Euler:}\qquad \begin{aligned} \mathbf{x}_{k+1}^{(p)} &= \mathbf{x}_k+h\,\mathbf{f}(t_k,\mathbf{x}_k,\mathbf{u}_k) \\ \mathbf{x}_{k+1} &= \mathbf{x}_k+\dfrac{h}{2}\Big(\mathbf{f}(t_k,\mathbf{x}_k,\mathbf{u}_k) +\mathbf{f}(t_{k+1},\mathbf{x}_{k+1}^{(p)},\mathbf{u}_{k+1})\Big) \end{aligned} \]

These formulas are the backbone of time-domain simulation in control engineering. Every software tool (including Simulink solvers and ODE libraries) implements the same core ideas, often with adaptive step-size and higher-order variants.

In practice, one estimates observed order using step halving. If \( E(h) \) is an error metric (for example, maximum absolute error against a reference solution), then

\[ p_{\text{obs}} \approx \dfrac{\log\!\left(E(h)/E(h/2)\right)}{\log 2} \]

For Euler, \( p_{\text{obs}}\approx 1 \); for improved Euler, \( p_{\text{obs}}\approx 2 \).

The following Python program solves the scalar IVP \( y'=-2y+\sin t,\; y(0)=1 \), compares Euler and improved Euler, estimates convergence order, and plots the trajectories. The exact solution is

\[ y(t)=\dfrac{2\sin t-\cos t}{5}+\dfrac{6}{5}e^{-2t}. \]

Chapter15_Lesson1.py

# Chapter15_Lesson1.py
# Numerical simulation of an IVP using Euler and Improved Euler (Heun) methods
# Example IVP: y' = -2 y + sin(t), y(0) = 1

import numpy as np
import matplotlib.pyplot as plt


def f(t, y):
    return -2.0 * y + np.sin(t)


def exact_solution(t):
    # y(t) = (2 sin t - cos t)/5 + (6/5) e^{-2 t}
    return (2.0 * np.sin(t) - np.cos(t)) / 5.0 + (6.0 / 5.0) * np.exp(-2.0 * t)


def euler(fhandle, t0, tf, y0, h):
    n = int(round((tf - t0) / h))
    t = np.linspace(t0, tf, n + 1)
    y = np.zeros(n + 1)
    y[0] = y0
    for k in range(n):
        y[k + 1] = y[k] + h * fhandle(t[k], y[k])
    return t, y


def improved_euler(fhandle, t0, tf, y0, h):
    """Heun / explicit trapezoid method."""
    n = int(round((tf - t0) / h))
    t = np.linspace(t0, tf, n + 1)
    y = np.zeros(n + 1)
    y[0] = y0
    for k in range(n):
        slope1 = fhandle(t[k], y[k])
        y_pred = y[k] + h * slope1
        slope2 = fhandle(t[k] + h, y_pred)
        y[k + 1] = y[k] + 0.5 * h * (slope1 + slope2)
    return t, y


def max_abs_error(t, y_num):
    return np.max(np.abs(y_num - exact_solution(t)))


def estimate_order(err_h, err_h2):
    return np.log(err_h / err_h2) / np.log(2.0)


def main():
    t0, tf, y0 = 0.0, 5.0, 1.0
    hs = [0.2, 0.1, 0.05, 0.025]

    print("Step-size study (max-norm error on [0, 5])")
    print("h        Euler error      Heun error")
    e_errors = []
    h_errors = []
    for h in hs:
        t_e, y_e = euler(f, t0, tf, y0, h)
        t_h, y_h = improved_euler(f, t0, tf, y0, h)
        ee = max_abs_error(t_e, y_e)
        he = max_abs_error(t_h, y_h)
        e_errors.append(ee)
        h_errors.append(he)
        print(f"{h:<8.3f} {ee:<14.6e} {he:<14.6e}")

    print("\nEstimated orders (using successive halving)")
    for i in range(len(hs) - 1):
        p_e = estimate_order(e_errors[i], e_errors[i + 1])
        p_h = estimate_order(h_errors[i], h_errors[i + 1])
        print(f"h={hs[i]:.3f} -> {hs[i+1]:.3f}: Euler p≈{p_e:.3f}, Heun p≈{p_h:.3f}")

    # Plot one representative step size
    h_plot = 0.1
    t_e, y_e = euler(f, t0, tf, y0, h_plot)
    t_h, y_h = improved_euler(f, t0, tf, y0, h_plot)
    t_dense = np.linspace(t0, tf, 800)
    y_dense = exact_solution(t_dense)

    plt.figure(figsize=(9, 5))
    plt.plot(t_dense, y_dense, label="Exact")
    plt.plot(t_e, y_e, "o-", markersize=3, label="Euler (h=0.1)")
    plt.plot(t_h, y_h, "s-", markersize=3, label="Improved Euler (h=0.1)")
    plt.xlabel("t")
    plt.ylabel("y(t)")
    plt.title("IVP solution: Euler vs Improved Euler")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

    # Save CSV for cross-language comparison
    out_data = np.column_stack([t_h, y_h, exact_solution(t_h), np.abs(y_h - exact_solution(t_h))])
    np.savetxt("Chapter15_Lesson1_output.csv", out_data, delimiter=",",
               header="t,heun,exact,abs_error", comments="")


if __name__ == "__main__":
    main()

      

Python libraries used here are basic numerical tools (numpy, matplotlib). In later lessons, you can compare these custom implementations with scipy.integrate.solve_ivp.

The C++ version is written from scratch using the standard library only. It prints error tables and exports a CSV trajectory for external plotting (e.g., Python, MATLAB, Excel, or gnuplot).

Chapter15_Lesson1.cpp

// Chapter15_Lesson1.cpp
// Euler and Improved Euler (Heun) for y' = -2 y + sin(t), y(0)=1

#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <vector>

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

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

void euler(double t0, double tf, double y0, double h,
           std::vector<double>& t, std::vector<double>& y) {
    int n = static_cast<int>(std::round((tf - t0) / h));
    t.resize(n + 1);
    y.resize(n + 1);
    t[0] = t0;
    y[0] = y0;
    for (int k = 0; k < n; ++k) {
        t[k + 1] = t[k] + h;
        y[k + 1] = y[k] + h * f(t[k], y[k]);
    }
}

void improved_euler(double t0, double tf, double y0, double h,
                    std::vector<double>& t, std::vector<double>& y) {
    int n = static_cast<int>(std::round((tf - t0) / h));
    t.resize(n + 1);
    y.resize(n + 1);
    t[0] = t0;
    y[0] = y0;
    for (int k = 0; k < n; ++k) {
        t[k + 1] = t[k] + h;
        double s1 = f(t[k], y[k]);
        double ypred = y[k] + h * s1;
        double s2 = f(t[k + 1], ypred);
        y[k + 1] = y[k] + 0.5 * h * (s1 + s2);
    }
}

double max_abs_error(const std::vector<double>& t, const std::vector<double>& y) {
    double emax = 0.0;
    for (std::size_t k = 0; k < t.size(); ++k) {
        double e = std::abs(y[k] - exact_solution(t[k]));
        if (e > emax) {
            emax = e;
        }
    }
    return emax;
}

double estimate_order(double err_h, double err_h2) {
    return std::log(err_h / err_h2) / std::log(2.0);
}

int main() {
    const double t0 = 0.0;
    const double tf = 5.0;
    const double y0 = 1.0;
    std::vector<double> hs = {0.2, 0.1, 0.05, 0.025};

    std::vector<double> eErr, hErr;

    std::cout << "Step-size study (max error on [0,5])\n";
    std::cout << std::left << std::setw(10) << "h"
              << std::setw(18) << "Euler"
              << std::setw(18) << "ImprovedEuler" << "\n";

    for (double h : hs) {
        std::vector<double> tE, yE, tH, yH;
        euler(t0, tf, y0, h, tE, yE);
        improved_euler(t0, tf, y0, h, tH, yH);
        double ee = max_abs_error(tE, yE);
        double he = max_abs_error(tH, yH);
        eErr.push_back(ee);
        hErr.push_back(he);

        std::cout << std::left << std::setw(10) << h
                  << std::setw(18) << ee
                  << std::setw(18) << he << "\n";
    }

    std::cout << "\nEstimated order (successive halving)\n";
    for (std::size_t i = 0; i + 1 < hs.size(); ++i) {
        std::cout << "h=" << hs[i] << " -> " << hs[i + 1]
                  << " : Euler p~" << estimate_order(eErr[i], eErr[i + 1])
                  << ", Heun p~" << estimate_order(hErr[i], hErr[i + 1]) << "\n";
    }

    // Save representative trajectory to CSV
    std::vector<double> t, y;
    improved_euler(t0, tf, y0, 0.1, t, y);
    std::ofstream out("Chapter15_Lesson1_cpp_output.csv");
    out << "t,heun,exact,abs_error\n";
    for (std::size_t k = 0; k < t.size(); ++k) {
        double ex = exact_solution(t[k]);
        out << t[k] << "," << y[k] << "," << ex << "," << std::abs(y[k] - ex) << "\n";
    }
    out.close();

    return 0;
}

      

The Java implementation mirrors the same algorithms and demonstrates language-independent numerical logic. It also writes a CSV file for trajectory comparison.

Chapter15_Lesson1.java

// Chapter15_Lesson1.java
// Euler and Improved Euler (Heun) for y' = -2 y + sin(t), y(0)=1

import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.Locale;

public class Chapter15_Lesson1 {

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

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

    static double[][] euler(double t0, double tf, double y0, double h) {
        int n = (int)Math.round((tf - t0) / h);
        double[] t = new double[n + 1];
        double[] y = new double[n + 1];
        t[0] = t0;
        y[0] = y0;
        for (int k = 0; k < n; k++) {
            t[k + 1] = t[k] + h;
            y[k + 1] = y[k] + h * f(t[k], y[k]);
        }
        return new double[][] {t, y};
    }

    static double[][] improvedEuler(double t0, double tf, double y0, double h) {
        int n = (int)Math.round((tf - t0) / h);
        double[] t = new double[n + 1];
        double[] y = new double[n + 1];
        t[0] = t0;
        y[0] = y0;
        for (int k = 0; k < n; k++) {
            t[k + 1] = t[k] + h;
            double s1 = f(t[k], y[k]);
            double yPred = y[k] + h * s1;
            double s2 = f(t[k + 1], yPred);
            y[k + 1] = y[k] + 0.5 * h * (s1 + s2);
        }
        return new double[][] {t, y};
    }

    static double maxAbsError(double[] t, double[] y) {
        double emax = 0.0;
        for (int k = 0; k < t.length; k++) {
            double e = Math.abs(y[k] - exact(t[k]));
            if (e > emax) emax = e;
        }
        return emax;
    }

    static double estimateOrder(double errH, double errH2) {
        return Math.log(errH / errH2) / Math.log(2.0);
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);

        double t0 = 0.0, tf = 5.0, y0 = 1.0;
        double[] hs = {0.2, 0.1, 0.05, 0.025};
        double[] eErr = new double[hs.length];
        double[] hErr = new double[hs.length];

        System.out.println("Step-size study (max error on [0,5])");
        System.out.printf("%-10s %-16s %-16s%n", "h", "Euler", "ImprovedEuler");

        for (int i = 0; i < hs.length; i++) {
            double h = hs[i];
            double[][] outE = euler(t0, tf, y0, h);
            double[][] outH = improvedEuler(t0, tf, y0, h);

            eErr[i] = maxAbsError(outE[0], outE[1]);
            hErr[i] = maxAbsError(outH[0], outH[1]);

            System.out.printf("%-10.3f %-16.8e %-16.8e%n", h, eErr[i], hErr[i]);
        }

        System.out.println("\nEstimated order (successive halving)");
        for (int i = 0; i < hs.length - 1; i++) {
            double pe = estimateOrder(eErr[i], eErr[i + 1]);
            double ph = estimateOrder(hErr[i], hErr[i + 1]);
            System.out.printf("h=%.3f -> %.3f : Euler p~%.4f, Heun p~%.4f%n", hs[i], hs[i+1], pe, ph);
        }

        // Export one trajectory to CSV
        double[][] out = improvedEuler(t0, tf, y0, 0.1);
        try (PrintWriter pw = new PrintWriter(new FileWriter("Chapter15_Lesson1_java_output.csv"))) {
            pw.println("t,heun,exact,abs_error");
            for (int k = 0; k < out[0].length; k++) {
                double t = out[0][k];
                double y = out[1][k];
                double ex = exact(t);
                pw.printf(Locale.US, "%.10f,%.10f,%.10f,%.10f%n", t, y, ex, Math.abs(y - ex));
            }
        }
    }
}

      

The MATLAB script implements Euler and improved Euler directly, estimates observed order, and plots the numerical and exact trajectories. This is useful pedagogically before using built-in solvers such as ode45.

Chapter15_Lesson1.m

% Chapter15_Lesson1.m
% Euler and Improved Euler (Heun) for y' = -2 y + sin(t), y(0)=1

clear; clc; close all;

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

t0 = 0; tf = 5; y0 = 1;
hs = [0.2, 0.1, 0.05, 0.025];

eErr = zeros(size(hs));
hErr = zeros(size(hs));

fprintf('Step-size study (max error on [0,5])\n');
fprintf('%-10s %-16s %-16s\n', 'h', 'Euler', 'ImprovedEuler');

for i = 1:numel(hs)
    h = hs(i);
    [tE, yE] = eulerMethod(f, t0, tf, y0, h);
    [tH, yH] = improvedEulerMethod(f, t0, tf, y0, h);

    eErr(i) = max(abs(yE - exact(tE)));
    hErr(i) = max(abs(yH - exact(tH)));

    fprintf('%-10.3f %-16.8e %-16.8e\n', h, eErr(i), hErr(i));
end

fprintf('\nEstimated order (successive halving)\n');
for i = 1:numel(hs)-1
    pE = log(eErr(i)/eErr(i+1))/log(2);
    pH = log(hErr(i)/hErr(i+1))/log(2);
    fprintf('h=%.3f -> %.3f : Euler p~%.4f, Heun p~%.4f\n', hs(i), hs(i+1), pE, pH);
end

hPlot = 0.1;
[tE, yE] = eulerMethod(f, t0, tf, y0, hPlot);
[tH, yH] = improvedEulerMethod(f, t0, tf, y0, hPlot);
tDense = linspace(t0, tf, 800);

figure;
plot(tDense, exact(tDense), 'LineWidth', 1.5); hold on;
plot(tE, yE, 'o-');
plot(tH, yH, 's-');
xlabel('t'); ylabel('y(t)');
title('Euler vs Improved Euler');
legend('Exact', 'Euler (h=0.1)', 'Improved Euler (h=0.1)', 'Location', 'best');
grid on;

T = table(tH(:), yH(:), exact(tH(:))', abs(yH(:)-exact(tH(:))'), ...
    'VariableNames', {'t', 'Heun', 'Exact', 'AbsError'});
writetable(T, 'Chapter15_Lesson1_matlab_output.csv');

function [t, y] = eulerMethod(f, t0, tf, y0, h)
    n = round((tf - t0)/h);
    t = linspace(t0, tf, n+1);
    y = zeros(1, n+1);
    y(1) = y0;
    for k = 1:n
        y(k+1) = y(k) + h * f(t(k), y(k));
    end
end

function [t, y] = improvedEulerMethod(f, t0, tf, y0, h)
    n = round((tf - t0)/h);
    t = linspace(t0, tf, n+1);
    y = zeros(1, n+1);
    y(1) = y0;
    for k = 1:n
        s1 = f(t(k), y(k));
        yPred = y(k) + h*s1;
        s2 = f(t(k+1), yPred);
        y(k+1) = y(k) + 0.5*h*(s1 + s2);
    end
end

      

Simulink note (conceptual setup): To replicate this IVP in Simulink, use a block diagram with a Clock block, a Trigonometric Function block for \( \sin t \), a gain of \( -2 \) acting on the state output, a summation block to form \( -2y+\sin t \), and an Integrator block with initial condition \( 1 \). In solver settings, you can compare a fixed-step Euler solver (if available) versus a higher-order solver and observe step-size effects directly.

The following Wolfram Language code implements the same methods and computes the experimental order of convergence. The file is named with the required convention and can be stored as a notebook file.

Chapter15_Lesson1.nb

(* Chapter15_Lesson1.nb *)
(* Plain-text Wolfram Language notebook-style code for Euler and Improved Euler (Heun) *)

ClearAll["Global`*"];

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

eulerMethod[f_, {t0_, tf_}, y0_, h_] := Module[
  {n, t, y, k},
  n = Round[(tf - t0)/h];
  t = Table[t0 + k h, {k, 0, n}];
  y = ConstantArray[0.0, n + 1];
  y[[1]] = y0;
  For[k = 1, k <= n, k++,
    y[[k + 1]] = y[[k]] + h f[t[[k]], y[[k]]];
  ];
  {t, y}
];

improvedEulerMethod[f_, {t0_, tf_}, y0_, h_] := Module[
  {n, t, y, k, s1, yPred, s2},
  n = Round[(tf - t0)/h];
  t = Table[t0 + k h, {k, 0, n}];
  y = ConstantArray[0.0, n + 1];
  y[[1]] = y0;
  For[k = 1, k <= n, k++,
    s1 = f[t[[k]], y[[k]]];
    yPred = y[[k]] + h s1;
    s2 = f[t[[k + 1]], yPred];
    y[[k + 1]] = y[[k]] + (h/2) (s1 + s2);
  ];
  {t, y}
];

maxAbsError[t_List, y_List] := Max[Abs[y - (exact /@ t)]];
estimateOrder[e1_, e2_] := Log[e1/e2]/Log[2];

t0 = 0; tf = 5; y0 = 1;
hs = {0.2, 0.1, 0.05, 0.025};

eErr = {}; hErr = {};
Print["Step-size study (max error on [0,5])"];
Do[
  {tE, yE} = eulerMethod[f, {t0, tf}, y0, h];
  {tH, yH} = improvedEulerMethod[f, {t0, tf}, y0, h];
  ee = maxAbsError[tE, yE];
  he = maxAbsError[tH, yH];
  AppendTo[eErr, ee];
  AppendTo[hErr, he];
  Print[Row[{"h=", NumberForm[h, {4, 3}], " | Euler=", ScientificForm[ee], " | Heun=", ScientificForm[he]}]];
, {h, hs}];

Print["Estimated order (successive halving)"];
Do[
  Print[
    Row[{
      "h=", hs[[i]], " -> ", hs[[i + 1]],
      " : Euler p~", NumberForm[estimateOrder[eErr[[i]], eErr[[i + 1]]], {4, 3}],
      ", Heun p~", NumberForm[estimateOrder[hErr[[i]], hErr[[i + 1]]], {4, 3}]
    }]
  ];
, {i, 1, Length[hs] - 1}];

{tEPlot, yEPlot} = eulerMethod[f, {t0, tf}, y0, 0.1];
{tHPlot, yHPlot} = improvedEulerMethod[f, {t0, tf}, y0, 0.1];

p1 = Plot[exact[t], {t, 0, 5}, PlotLegends -> {"Exact"}];
p2 = ListLinePlot[
  {
    Transpose[{tEPlot, yEPlot}],
    Transpose[{tHPlot, yHPlot}]
  },
  PlotLegends -> {"Euler (h=0.1)", "Improved Euler (h=0.1)"}
];
Show[p1, p2, GridLines -> Automatic, AxesLabel -> {"t", "y"}]

      

Problem 1 (Euler local truncation error): For the scalar IVP \( y'=f(t,y) \) with sufficiently smooth exact solution \( y(t) \), prove that the Euler local truncation error satisfies \( \eta_{k+1}^{(E)}=O(h^2) \).

Solution: Using Taylor expansion about \( t_k \),

\[ y(t_{k+1}) = y(t_k)+h y'(t_k)+\dfrac{h^2}{2}y''(\xi_k) \]

and since \( y'(t_k)=f(t_k,y(t_k)) \),

\[ \eta_{k+1}^{(E)} = y(t_{k+1})-y(t_k)-h f(t_k,y(t_k)) = \dfrac{h^2}{2}y''(\xi_k). \]

If \( |y''(t)| \le M \) on the interval, then \( |\eta_{k+1}^{(E)}| \le \dfrac{M}{2}h^2 \), hence \( O(h^2) \).

Problem 2 (Global error recursion): Assume \( f \) is Lipschitz in \( y \) with constant \( L \), and Euler LTE satisfies \( |\eta_{k+1}^{(E)}| \le Ch^2 \). Show the Euler global error is \( O(h) \).

Solution: Let \( e_k=y(t_k)-y_k \). Then

\[ e_{k+1} = e_k + h\left(f(t_k,y(t_k))-f(t_k,y_k)\right)+\eta_{k+1}^{(E)} \]

so

\[ |e_{k+1}| \le (1+hL)|e_k| + Ch^2. \]

With \( e_0=0 \), discrete Grönwall gives

\[ |e_k| \le \dfrac{C}{L}\left(e^{L(t_k-t_0)}-1\right)h. \]

Therefore \( \max_k |e_k| = O(h) \).

Problem 3 (Improved Euler one-step computation): Apply one step of improved Euler to \( y'=t-y \), \( y(0)=1 \), with \( h=0.1 \). Also compare with Euler.

Solution: Euler slope at \( (0,1) \) is

\[ s_1 = f(0,1)=0-1=-1 \]

Euler step:

\[ y_1^{(E)} = 1 + 0.1(-1)=0.9 \]

Improved Euler predictor:

\[ y_1^{(p)} = 1+0.1(-1)=0.9 \]

End slope using the predicted value:

\[ s_2=f(0.1,0.9)=0.1-0.9=-0.8 \]

Corrected step:

\[ y_1 = 1+\dfrac{0.1}{2}(s_1+s_2)=1+0.05(-1.8)=0.91. \]

The exact solution is \( y(t)=t-1+2e^{-t} \), so \( y(0.1)\approx 0.909674 \). Improved Euler is closer than Euler in this step.

Problem 4 (Stability-limited step size): For the decay system \( y'=-10y \), determine a stable fixed step range for Euler and improved Euler on the negative real axis.

Solution: Here \( \lambda=-10 \), so \( z=h\lambda=-10h \). For both Euler and improved Euler on the negative real axis, the stability interval is

\[ 0 < -h\lambda < 2 \]

hence

\[ 0 < 10h < 2 \quad \Rightarrow \quad 0 < h < 0.2. \]

So any fixed step \( h \ge 0.2 \) is at or beyond the stability boundary and is not acceptable for robust simulation.

Problem 5 (State update for a mass–spring–damper system): Consider \( m\ddot{q}+c\dot{q}+kq=0 \). Let \( x_1=q \), \( x_2=\dot{q} \). Derive the Euler updates.

Solution: The first-order model is

\[ \dot{x}_1=x_2, \qquad \dot{x}_2=-\dfrac{k}{m}x_1-\dfrac{c}{m}x_2. \]

Euler updates at \( t_k \) are

\[ x_{1,k+1}=x_{1,k}+h x_{2,k} \]

\[ x_{2,k+1}=x_{2,k}+h\left(-\dfrac{k}{m}x_{1,k}-\dfrac{c}{m}x_{2,k}\right). \]

These formulas directly simulate displacement and velocity. Improved Euler uses the same model but computes a predictor state first and then averages slopes.

We formulated IVPs for dynamic systems, derived Euler and improved Euler from Taylor ideas, and proved their accuracy orders using local truncation error and global error recursion. We also analyzed numerical stability using the test equation \( y'=\lambda y \), showing how the step size must respect the dynamics of the system. The provided implementations in Python, C++, Java, MATLAB/Simulink, and Mathematica establish a reusable template for later lessons on Runge–Kutta methods, stiffness, and adaptive time stepping.