Numerical Integration and Differential Equations

MATLAB's ODE solvers integrate dy/dt = f(t, y) from initial conditions:

Concept

Integrate sampled data (trapezoidal rule) cumtrapz(x, y)

Running integral (e.g. velocity from acceleration) integral(@f, a, b)

Adaptive integration of a function (infinite limits allowed) integral2, integral3 Double and triple integrals Ordinary Differential Equations

MATLAB's ODE solvers integrate dy/dt = f(t, y) from initial conditions:

Default: explicit Runge–Kutta 4(5), non-stiff problems Lower accuracy, crude tolerances ode15s, ode23s

Stiff problems (very different time scales: chemical kinetics, circuits with tiny time constants)

Tolerances (RelTol, AbsTol), events, max step Higher-order ODEs are rewritten as first-order systems: x'' + 2ζωₙx' + ωₙ²x = 0 becomes s₁' = s₂, s₂' = −2ζωₙs₂ − ωₙ²s₁.

The solvers choose adaptive time steps; the returned t is not equally spaced (pass a vector of times as tspan to get output at those times).

Fixed-step methods like Euler (yₙ₊₁ = yₙ + h·f) are simple and used in embedded controllers, but their error shrinks only in proportion to h.

[t, y] = ode45(@(t, y) f(t, y), [t0 tf], y0, options);

Why It Matters

Dynamics of vehicles, circuits, thermal systems, batteries and chemical reactors are ODEs. Simulating them before building hardware is the heart of model-based design (Simulink solves the same equations graphically).

Syntax

area = trapz(t, v);
[t, y] = ode45(@(t, y) -y / tau + u / tau, [0 5], 0);
opts = odeset('RelTol', 1e-6, 'Events', @stopWhenEmpty);

Example

23.1 Integration, an RC circuit, a damped oscillator and Euler's method

% numerical integration
x = linspace(0, pi, 101);
fprintf('trapz(sin) = %.6f, integral(sin) = %.10f (exact 2)\n', trapz(x, sin(x)), integral(@sin, 0, pi))
fprintf('integral of exp(-x^2) over the real line = %.8f (sqrt(pi) = %.8f)\n', integral(@(x) exp(-x .^ 2), -Inf, Inf), sqrt(pi))

% RC circuit charging: dv/dt = (Vin - v) / (R*C)
Rr = 1e3; C = 1e-3; Vin = 5;
rc = @(t, v) (Vin - v) / (Rr * C);
[t, v] = ode45(rc, [0 5], 0);
fprintf('RC: v(1 s) = %.4f V (exact %.4f), v(5 s) = %.4f V\n', interp1(t, v, 1), Vin * (1 - exp(-1)), v(end))

% mass-spring-damper as a first-order system: x'' + 2*zeta*wn x' + wn^2 x = 0
wn = 2 * pi; zeta = 0.1;
f = @(t, s) [s(2); -2 * zeta * wn * s(2) - wn ^ 2 * s(1)];
opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[t2, s2] = ode45(f, 0:0.001:4, [1; 0], opts);   % a time vector fixes the output grid
pk = find(s2(2:end-1, 1) > s2(1:end-2, 1) & s2(2:end-1, 1) > s2(3:end, 1)) + 1;
fprintf('damped oscillation: first peaks at t = %s s, amplitudes %s\n', mat2str(t2(pk(1:3))', 3), mat2str(s2(pk(1:3), 1)', 3))
fprintf('logarithmic decrement %.4f (theory 2*pi*zeta/sqrt(1-zeta^2) = %.4f)\n', ...
        log(s2(pk(1), 1) / s2(pk(2), 1)), 2 * pi * zeta / sqrt(1 - zeta ^
…

Expected Result

integral handled an infinite interval (the Gaussian integral √π).

The RC voltage after one time constant (RC = 1 s) agreed with 5·(1 − e⁻¹) V to within 0.0003 V — the gap comes from ode45's default tolerances (tighten them with odeset('RelTol', ...)).

The damped oscillator's successive peaks shrink by the logarithmic decrement predicted by theory.

Euler's error dropped roughly tenfold for each tenfold smaller step — first-order accuracy.

trapz(sin) = 1.999836, integral(sin) = 2.0000000000 (exact 2)
integral of exp(-x^2) over the real line = 1.77245385 (sqrt(pi) = 1.77245385)
RC: v(1 s) = 3.1603 V (exact 3.1606), v(5 s) = 4.9663 V
damped oscillation: first peaks at t = [1.01 2.01 3.02] s, amplitudes [0.532 0.283 0.15]
logarithmic decrement 0.6315 (theory 2*pi*zeta/sqrt(1-zeta^2) = 0.6315)
Euler h = 0.5 : v(1) = 3.7500 (error 0.5894)
Euler h = 0.1 : v(1) = 3.2566 (error 0.0960)
Euler h = 0.01: v(1) = 3.1698 (error 0.0092)

Common Mistakes

Writing the ODE function for a second-order system without converting it to first order.

Using ode45 on stiff problems (extremely slow) — switch to ode15s.

Assuming the solver's t is equally spaced.

Loose default tolerances for sensitive problems.

Best Practices

Check the solution against an analytic case or a smaller tolerance.

Pass parameters by capturing them in an anonymous function.

Use Events to stop at a condition (battery empty, object hits the ground).

Interview Relevance

Medium–high for engineering: "how do you solve an ODE in MATLAB?", "what is a stiff problem?", "how do you convert a second-order ODE to first order?", "ode45 vs Euler?".

Practice Questions

Integrate an acceleration log twice with cumtrapz to get position.

Simulate a falling object with air drag until it hits the ground using an event function.

Simulate an RLC circuit and find its resonant frequency from the response.

Compare ode45 and ode15s run times on the Van der Pol equation with μ = 1000.

Quick Quiz

Q1. (Multiple choice) Which solver is designed for stiff ODEs? A) ode45 B) ode23 C) ode15s D) trapz Answer: C.

Q2. (Predict the output) RC charging with Vin = 5 V, RC = 1 s: v(1 s)?

Answer: ≈ 3.1606 V (verified).

Q3. (True/False) trapz needs a function handle.

Answer: False. It integrates sampled data; integral takes a handle.

Q4. (Find the bug) ode45(@(t, x) [x(2); -k*x(1)], [0 10], 1) errors.

Answer: Two state variables need two initial conditions: [1; 0].

Coding Challenge

Simulate the logistic population model dP/dt = 0.5·P·(1 − P/1000) from P(0) = 10 for 30 time units with ode45, and print P at t = 10, 20, 30.

f = @(t, P) 0.5 * P .* (1 - P / 1000);
[t, P] = ode45(f, [0 10 20 30], 10);
fprintf('P(%d) = %.1f\n', [t(2:end)'; P(2:end)'])