# Solving a Mass-Spring IVP in MATLAB
This page gets you from the equation
$m\,y''+b\,y'+k\,y=f(t),\qquad y(t_0)=y_0,\quad y'(t_0)=v_0$
to a plot of its solution, using MATLAB's built-in solver `ode45`. You will not write any code. You will open one script, change a few numbers at the top, and press Run. The rest of the page explains what the script and the solver are doing, how to check that an answer can be trusted, and what each standard case should look like so you can compare your output against a reference.
You need MATLAB (any recent version; `ode45` is built in and needs no toolbox) and the four scripts collected under [[#Full codes]] at the bottom of this page. Save each one as a `.m` file with the name shown, all in the same folder.
> [!info]- What you will have at the end
> - `simulateSpring.m` solves and plots any mass-spring problem you describe in its first section.
> - `checkSolver.m` shows two checks for deciding whether a numerical answer is right.
> - `underTheHood.m` shows how a time stepper works and where `ode45` puts its steps.
> - `missedPush.m` shows the most common way a solver can quietly give a wrong answer, and the fix.
# The idea
## What "solving numerically" means
For the cases you met in a differential equations course, $y(t)$ comes out as a formula. Once the forcing, the coefficients, or the model get more realistic, that stops being possible. A numerical solver gives you something different. It gives you a list of times $t_0<t_1<t_2<\dots$ and, at each one, an approximate position $y_n\approx y(t_n)$ and velocity $v_n\approx y'(t_n)$. Plotting that list is plotting the solution.
The solver builds the list one step at a time. It knows the state now, it can compute how fast the state is changing, and it uses that to move a short time $h$ forward. Then it repeats.
## Rewrite the equation as a first-order system
Solvers only take one kind of input, a first-order system $\vec Y\,'=\vec F(t,\vec Y)$. A second-order equation is turned into one by giving the velocity its own name,
$v=y'.$
Then $y''=v'$, and solving the original equation for $y''$ gives
$\begin{aligned}
y' &= v,\\
v' &= \frac{f(t)-b\,v-k\,y}{m}.
\end{aligned}$
Stack the two unknowns into a **state vector** $\vec Y=\begin{bmatrix}y\\ v\end{bmatrix}$. The system says how fast each entry of the state changes, given the time and the current state,
$\vec Y\,'=\vec F(t,\vec Y)=\begin{bmatrix} v\\[2pt] \dfrac{f(t)-b\,v-k\,y}{m}\end{bmatrix},\qquad \vec Y(t_0)=\begin{bmatrix}y_0\\ v_0\end{bmatrix}.$
That is all a solver needs. Two initial conditions become one initial state, and one second-order equation becomes two first-order equations. Any $n$th-order equation can be handled the same way by carrying $y,\,y',\,\dots,\,y^{(n-1)}$ in the state.
# The script, piece by piece
`simulateSpring.m` has five numbered sections. Only Section 1 is meant to be edited. The full listing is under [[#Full codes]].
## Section 1. The problem data
> [!example]- MATLAB - Section 1, problem data (edit these)
> *Every number that defines the problem lives here. Change these and nothing else.*
> ```matlab
> m = 1; % mass
> b = 0.5; % damping coefficient
> k = 1; % spring constant
> f = @(t) 0; % external forcing f(t); @(t) 0 means no forcing
>
> y0 = 1; % initial position y(t0)
> v0 = 0; % initial velocity y'(t0)
> tspan = [0 30]; % [start time, end time]
> ```
The forcing is written as an **anonymous function**, MATLAB's way of writing a short formula in $t$. The `@(t)` says "a function of $t$," and the expression after it is the formula. Some forcings you might want are below.
| Forcing | How to write it |
|---|---|
| none | `f = @(t) 0;` |
| constant push of size 2 | `f = @(t) 2;` |
| $F_0\cos(\omega t)$ with $F_0=3$, $\omega=0.8$ | `f = @(t) 3*cos(0.8*t);` |
| push of size 5 between $t=10$ and $t=10.5$ | `f = @(t) 5*(t >= 10 & t <= 10.5);` |
| decaying kick $e^{-t}\sin(2t)$ | `f = @(t) exp(-t).*sin(2*t);` |
The comparison `(t >= 10 & t <= 10.5)` is $1$ when true and $0$ when false, which is how an on/off forcing is written in one line. If you use a forcing with a short burst like this, read [[#When things go wrong]] first.
## Section 2. The first-order system
> [!example]- MATLAB - Section 2, the right-hand side $\vec F(t,\vec Y)$
> *This is the system from [[#Rewrite the equation as a first-order system]], typed in. `Y(1)` is $y$ and `Y(2)` is $v$.*
> ```matlab
> rhs = @(t, Y) [ Y(2);
> ( f(t) - b*Y(2) - k*Y(1) ) / m ];
> ```
`rhs` is a function of two inputs, the time `t` and the state `Y`, and it returns a column of two numbers, $y'$ and $v'$. The semicolon inside the brackets stacks the two entries vertically, and `ode45` requires the result to be a column. Because `m`, `b`, `k`, and `f` already exist when this line runs, `rhs` remembers their values.
## Section 3. Solver settings
> [!example]- MATLAB - Section 3, tolerances
> *How much error `ode45` is allowed to commit on each step.*
> ```matlab
> opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
> ```
`RelTol` is the allowed error relative to the size of the solution, so `1e-8` asks for roughly eight correct digits. `AbsTol` is a floor that matters when the solution is near zero, where "relative to the size" stops meaning anything. If you leave the options out, MATLAB uses `RelTol = 1e-3` and `AbsTol = 1e-6`. Those defaults are loose, and [[#Can you trust the answer]] shows what they can do over a long run. For a two-equation system, tight tolerances cost a fraction of a second, so the script sets them tight.
## Section 4. Solve
> [!example]- MATLAB - Section 4, the call to `ode45`
> *Hand over the system, the time interval, the initial state as a column, and the settings.*
> ```matlab
> [t, Y] = ode45(rhs, tspan, [y0; v0], opts);
> y = Y(:, 1); % position at each output time
> v = Y(:, 2); % velocity at each output time
> ```
`ode45` returns two things. `t` is a column of times. `Y` has one row per time and one column per state variable, so `Y(:,1)` is the position and `Y(:,2)` is the velocity.
Two things about `t` surprise people. The times are **not evenly spaced**, because `ode45` chooses its own step sizes (see [[#Adaptive steps, which is what ode45 adds]]). And the number of output times is **not** the number of steps taken, because MATLAB fills in a few extra points inside each step so plots look smooth. If you want the solution at particular times, give `tspan` as a list instead of two endpoints, for example `tspan = linspace(0, 30, 301);` for every $0.1$. The solver still picks its own steps and then reports the solution at the times you asked for.
## Section 5. Plot
> [!example]- MATLAB - Section 5, position and phase plane
> *Left panel is $y(t)$. Right panel is $v$ against $y$, with a dot at the starting state.*
> ```matlab
> figure('Position', [100 100 1000 400]);
>
> subplot(1, 2, 1);
> plot(t, y, 'LineWidth', 1.5);
> grid on; xlim([tspan(1) tspan(end)]);
> xlabel('t'); ylabel('y(t)');
> title(sprintf('m = %g, b = %g, k = %g, f = %s', m, b, k, func2str(f)));
>
> subplot(1, 2, 2);
> plot(y, v, 'LineWidth', 1.5); hold on;
> plot(y(1), v(1), 'ko', 'MarkerFaceColor', 'k'); % starting point
> hold off; grid on;
> xlabel('y'); ylabel('v = y''');
> title('Phase plane (dot = start)');
> ```
The **phase plane** plots the state itself, velocity against position, with time left implicit. It shows the type of motion at a glance. A closed loop means the motion repeats forever, an inward spiral means oscillation that dies out, and a curve that falls straight into the origin means the mass returns to rest without swinging.
## What you should see
With the values in Section 1 as given ($m=1$, $b=0.5$, $k=1$, no forcing, released from $y=1$ at rest), pressing Run prints
```text
ode45 returned 1345 output times on [0, 30].
Final state: y = -0.000495616, v = 0.000398894
```
and draws
![[ivp-tutorial_underdamped.png]]
The final state should match to at least four digits. The number of output times may differ between MATLAB versions, and it does not matter. All figures on this page were made by running these exact scripts in GNU Octave, which runs the same MATLAB code, so your fonts, colors, and the way the forcing is printed in the title will look slightly different. The curves and the printed numbers are what to compare.
# Under the hood
You do not need this section to use the script. It explains what happens between pressing Run and getting the list of times and states, which is what you need to read a solver's warnings and to trust its output.
## A time stepper
Every method here does the same thing. Start from the known state $\vec Y_n$ at time $t_n$, use $\vec F$ to estimate how the state changes over a short time $h$, and land on $\vec Y_{n+1}$ at $t_{n+1}=t_n+h$.
The simplest choice is **Euler's method**. Compute the slope at the start of the step and follow it in a straight line,
$\vec Y_{n+1}=\vec Y_n+h\,\vec F(t_n,\vec Y_n).$
This is the definition of the derivative read backwards. Its weakness is that the slope changes during the step and Euler ignores that. The error made over a whole run shrinks only in proportion to $h$, so halving the step only halves the error.
## Better slopes with RK4
The classical **fourth-order Runge-Kutta method (RK4)** samples the slope four times across the step and takes a weighted average.
$\begin{aligned}
\vec k_1 &= \vec F(t_n,\ \vec Y_n)\\
\vec k_2 &= \vec F\!\left(t_n+\tfrac h2,\ \vec Y_n+\tfrac h2\vec k_1\right)\\
\vec k_3 &= \vec F\!\left(t_n+\tfrac h2,\ \vec Y_n+\tfrac h2\vec k_2\right)\\
\vec k_4 &= \vec F\!\left(t_n+h,\ \vec Y_n+h\,\vec k_3\right)\\[4pt]
\vec Y_{n+1} &= \vec Y_n+\frac h6\left(\vec k_1+2\vec k_2+2\vec k_3+\vec k_4\right)
\end{aligned}$
The first slope is Euler's. The second and third look at the midpoint of the step, and the fourth looks at the end. The weights $1,2,2,1$ are chosen so that the step agrees with the exact solution's Taylor series through the $h^4$ term. Over a whole run the error shrinks like $h^4$, so halving the step cuts the error by about $16$.
> [!example]- MATLAB - one Euler step and one RK4 step (reading only)
> *Euler, one slope.*
> ```matlab
> YE(:,n+1) = YE(:,n) + h * rhs(t(n), YE(:,n));
> ```
> *RK4, four slopes and their weighted average.*
> ```matlab
> k1 = rhs(t(n), YR(:,n));
> k2 = rhs(t(n) + h/2, YR(:,n) + h/2*k1);
> k3 = rhs(t(n) + h/2, YR(:,n) + h/2*k2);
> k4 = rhs(t(n) + h, YR(:,n) + h*k3);
> YR(:,n+1) = YR(:,n) + h/6 * (k1 + 2*k2 + 2*k3 + k4);
> ```
The left panel below is Part A of `underTheHood.m`. Both methods take the same $80$ steps of size $h=0.25$ on the underdamped spring.
```text
Part A, step h = 0.25 (80 steps)
Euler max error: 2.58e-01
RK4 max error: 5.21e-05
```
![[ivp-tutorial_underTheHood.png]]
Euler's oscillation dies out far too slowly and drifts out of phase. RK4 sits on the exact curve, about $5000$ times more accurate, for four times the work per step.
## Adaptive steps, which is what ode45 adds
A fixed step has to be small enough for the hardest part of the run, which wastes effort everywhere else, and you have to guess it in advance. `ode45` picks the step for you as it goes.
On every step it computes **two** estimates of $\vec Y_{n+1}$ from the same set of slope samples, one accurate to fourth order and one to fifth order. (The method is the Dormand-Prince pair, a close relative of RK4 that uses six new slope samples per step.) The two estimates agree closely when the step is small enough and disagree when it is not, so their difference $\vec e$ is a free estimate of the error just committed. The step is accepted only if every component passes
$|e_i|\le\max\big(\texttt{RelTol}\cdot|y_i|,\ \texttt{AbsTol}\big).$
If it fails, the step is thrown away and retried with a smaller $h$. If it passes easily, the next step is allowed to be larger. The number in `ode45` refers to this 4/5 pair.
The right panels above are Part B of `underTheHood.m`. A lightly damped spring is released from $y=1$ and then pushed with $f=5$ for $10\le t\le10.5$. Each dot is one step `ode45` actually took. During the smooth ringing the steps stay near $0.2$. Where the push switches on and off, the forcing jumps, and the steps drop by four orders of magnitude to get through the corner, then grow back right away.
```text
Part B: 208 steps, smallest 2.45e-05, largest 2.17e-01
```
# Can you trust the answer
A solver always returns numbers, whether or not they are right. Two checks tell you whether they are, and neither one needs new code. `checkSolver.m` runs both.
**Check A. Compare against a case with a known exact answer.** Before trusting the solver on a problem you cannot solve by hand, run it on one you can. The resonance problem $y''+y=\cos t$ with $y(0)=y'(0)=0$ has exact solution $y=\tfrac t2\sin t$. The left panel shows the error over time for three tolerances.
```text
Check A: resonance, compared with y = (t/2) sin t
RelTol output times max |error|
1e-03 325 8.48e-02
1e-06 1173 4.20e-05
1e-09 4721 3.61e-08
```
![[ivp-tutorial_checkSolver.png]]
The error tracks the tolerance. Each factor of $1000$ in `RelTol` buys about three more correct digits, at the cost of roughly four times as many output times.
**Check B. Tighten the tolerance and see whether the answer moves.** This is the check you can always do, even when no exact answer exists. Run the problem, then run it again with `RelTol` and `AbsTol` both $100$ times smaller. If the curves you care about do not visibly change, the first answer was good enough. If they do change, keep tightening until they stop.
The right panel shows why this matters. An undamped spring $y''+y=0$ released from $y=1$ should circle the origin of the phase plane at distance $1$ forever. With MATLAB's default tolerances it loses about $3\%$ of its amplitude by $t=1000$, which looks exactly like a small amount of damping that is not in the equation. With tight tolerances it stays at $1$ to six digits.
```text
Check B: undamped spring on [0, 1000], distance from origin at the end
default tolerances: 0.9701 (should be 1)
tight tolerances: 0.99999938
```
Your MATLAB numbers for the default run may differ a little from $0.9701$. The point is that they are visibly not $1$.
# Reference cases
Each case below uses $m=1$ and $k=1$, so the natural frequency is $\omega_0=\sqrt{k/m}=1$ and the period of free oscillation is $2\pi$. Damping divides into regimes at $b=2\sqrt{mk}=2$. Set Section 1 of `simulateSpring.m` to the values in a row, press Run, and compare your printed final state and your figure to the ones here.
| Case | $b$ | `f` | $y_0$ | $v_0$ | `tspan` | Final $y$ | Final $v$ |
|---|---|---|---|---|---|---|---|
| Undamped | `0` | `@(t) 0` | `1` | `0` | `[0 30]` | $0.1543$ | $0.9880$ |
| Underdamped | `0.5` | `@(t) 0` | `1` | `0` | `[0 30]` | $-4.956\times10^{-4}$ | $3.989\times10^{-4}$ |
| Critically damped | `2` | `@(t) 0` | `1` | `0` | `[0 20]` | $4.329\times10^{-8}$ | $-4.123\times10^{-8}$ |
| Overdamped | `4` | `@(t) 0` | `1` | `0` | `[0 20]` | $5.070\times10^{-3}$ | $-1.358\times10^{-3}$ |
| Driven, damped | `0.5` | `@(t) cos(0.8*t)` | `0` | `0` | `[0 60]` | $-1.857$ | $0.05667$ |
| Resonance | `0` | `@(t) cos(t)` | `0` | `0` | `[0 50]` | $-6.559$ | $23.99$ |
| Beating | `0` | `@(t) cos(0.9*t)` | `0` | `0` | `[0 150]` | $-8.923$ | $-4.181$ |
Where an exact solution is easy to write down, the table values agree with it. For example, the undamped case is $y=\cos t$, so $y(30)=\cos 30=0.1543$, and the resonance case is $y=\tfrac t2\sin t$, so $y(50)=25\sin 50=-6.559$.
## Undamped, $b=0$
![[ivp-tutorial_undamped.png]]
The solution is $y=\cos t$. With no damping the energy is conserved, so the phase plane is a closed circle traced clockwise over and over.
## Underdamped, $0<b<2$
![[ivp-tutorial_underdamped.png]]
The mass still oscillates, but each swing is smaller than the last. The amplitude decays like $e^{-bt/(2m)}$ and the phase plane spirals into the origin.
## Critically damped, $b=2$
![[ivp-tutorial_critical.png]]
The fastest return to rest without overshooting. The solution is $y=(1+t)e^{-t}$. In the phase plane the curve drops straight into the origin without circling it.
## Overdamped, $b>2$
![[ivp-tutorial_overdamped.png]]
Again no oscillation, but the return is slower than critical damping because the heavy damping resists the motion. Compare the two $y(t)$ panels. At $t=10$ the critically damped mass is essentially at rest while the overdamped one is still visibly creeping in.
## Driven and damped, $f=\cos(0.8t)$
![[ivp-tutorial_steadystate.png]]
Starting from rest, the response is a **transient** that dies out plus a **steady-state** oscillation at the driving frequency $0.8$ that does not. After about $t=20$ only the steady state is left, with amplitude
$\frac{1}{\sqrt{(k-m\omega^2)^2+(b\,\omega)^2}}=\frac{1}{\sqrt{0.36^2+0.4^2}}\approx1.858.$
In the phase plane the spiral settles onto a closed loop, the steady-state orbit.
## Resonance, $b=0$ and $f=\cos t$
![[ivp-tutorial_resonance.png]]
Driving an undamped spring exactly at its natural frequency adds energy on every swing. The exact solution is $y=\tfrac t2\sin t$, so the amplitude grows linearly without bound and the phase plane spirals outward.
## Beating, $b=0$ and $f=\cos(0.9t)$
![[ivp-tutorial_beating.png]]
Driving close to, but not at, the natural frequency. The exact solution is
$y=\frac{\cos(0.9t)-\cos t}{1-0.81}=\frac{2}{0.19}\,\sin(0.95t)\,\sin(0.05t),$
a fast oscillation inside a slow envelope. The amplitude swells to about $10.5$ and dies back to zero every $\pi/0.05\approx62.8$ time units. The phase plane fills in a disk instead of tracing one curve.
# When things go wrong
## A short forcing gets skipped
This is the one failure that gives no warning. Start a spring at rest and push it with $f=5$ for $10\le t\le10.5$. Before $t=10$ the solution is exactly zero, so `ode45`'s error estimate is zero and it grows its step as large as it is allowed. By default that limit is a tenth of the whole interval, here $3$. A step of size $3$ can begin before $t=10$ and end after $t=10.5$ without ever sampling $f$ inside the push, and the solver reports that nothing happened.
```text
No MaxStep: largest |y| = 0.0000
MaxStep = 0.1: largest |y| = 2.1341
```
![[ivp-tutorial_missedPush.png]]
The fix is to cap the step below the width of the shortest feature in the forcing, with the `MaxStep` option.
> [!example]- MATLAB - capping the step size
> *Add `MaxStep` to the options in Section 3. Here the push lasts $0.5$, so $0.1$ is safe.*
> ```matlab
> opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10, 'MaxStep', 0.1);
> ```
Whether MATLAB's `ode45` actually skips this particular push depends on its exact step sequence, so you may or may not reproduce the flat line. The risk is real either way, and `MaxStep` costs very little. Any time `f` switches on and off, set `MaxStep` smaller than the shortest on or off period.
## Other things to look for
- **An error saying the function must return a column vector**, or a similar message about sizes. The right-hand side is returning a row. Check that the two entries inside the brackets in Section 2 are separated by a semicolon, and that the initial state is written `[y0; v0]`.
- **An error inside `f`.** The forcing must be a function of `t` alone, written `@(t) ...`. A forcing written as a number, `f = 2`, will fail because Section 2 calls it as `f(t)`.
- **The answer changes when you change the tolerance.** The solution is not converged. Keep tightening until it stops changing (Check B).
- **`ode45` runs for a very long time on a simple-looking problem.** This usually means the problem is **stiff**, meaning it has two time scales that are very far apart. For example, $m=10^{-4}$, $b=1$, $k=1$ has one motion that dies out in about $10^{-4}$ time units and another that takes about one. The solution after the first instant is a slow, smooth decay, yet `ode45` needs tens of thousands of tiny steps to stay stable. MATLAB's `ode15s` is built for stiff problems and takes the same inputs, so replace `ode45` with `ode15s` in Section 4. For ordinary values of $m$, $b$, and $k$ this does not come up.
- **The plot looks wrong but there is no error.** Check the physics before the numerics. Is the sign of the damping right? Is the forcing divided by $m$? (It is, in Section 2.) Is `tspan` long enough to see what you are looking for?
# Full codes
Save each listing as a `.m` file with the name in its title, all in the same folder. Run a script by opening it and pressing Run, or by typing its name at the command prompt.
> [!example]- MATLAB - `simulateSpring.m` (the simulator; edit Section 1 only)
> ```matlab
> %% simulateSpring.m
> % Approximately solves the mass-spring initial value problem
> %
> % m y'' + b y' + k y = f(t), y(t0) = y0, y'(t0) = v0
> %
> % with MATLAB's built-in adaptive solver ode45. Edit the values in
> % Section 1, then press Run. Nothing below Section 1 needs to change.
>
> clear; close all;
>
> %% 1. Problem data (edit these)
> m = 1; % mass
> b = 0.5; % damping coefficient
> k = 1; % spring constant
> f = @(t) 0; % external forcing f(t); @(t) 0 means no forcing
>
> y0 = 1; % initial position y(t0)
> v0 = 0; % initial velocity y'(t0)
> tspan = [0 30]; % [start time, end time]
>
> %% 2. Rewrite as a first-order system
> % The state is Y = [y; v] with v = y'. Then
> % y' = v
> % v' = ( f(t) - b*v - k*y ) / m
> % Y(1) is position and Y(2) is velocity.
> rhs = @(t, Y) [ Y(2);
> ( f(t) - b*Y(2) - k*Y(1) ) / m ];
>
> %% 3. Solver settings
> % RelTol and AbsTol set how much error ode45 allows in each step.
> % (MATLAB's defaults are RelTol = 1e-3 and AbsTol = 1e-6, which are loose.)
> opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
>
> %% 4. Solve
> [t, Y] = ode45(rhs, tspan, [y0; v0], opts);
> y = Y(:, 1); % position at each output time
> v = Y(:, 2); % velocity at each output time
>
> fprintf('ode45 returned %d output times on [%g, %g].\n', numel(t), tspan(1), tspan(end));
> fprintf('Final state: y = %.6g, v = %.6g\n', y(end), v(end));
>
> %% 5. Plot
> figure('Position', [100 100 1000 400]);
>
> subplot(1, 2, 1);
> plot(t, y, 'LineWidth', 1.5);
> grid on; xlim([tspan(1) tspan(end)]);
> xlabel('t'); ylabel('y(t)');
> title(sprintf('m = %g, b = %g, k = %g, f = %s', m, b, k, func2str(f)));
>
> subplot(1, 2, 2);
> plot(y, v, 'LineWidth', 1.5); hold on;
> plot(y(1), v(1), 'ko', 'MarkerFaceColor', 'k'); % starting point
> hold off; grid on;
> xlabel('y'); ylabel('v = y''');
> title('Phase plane (dot = start)');
> ```
> [!example]- MATLAB - `checkSolver.m` (the two trust checks)
> ```matlab
> %% checkSolver.m
> % Two ways to decide whether an ode45 answer can be trusted.
> % Check A. Compare against a case whose exact solution is known.
> % Check B. Change the tolerance and see whether the answer moves.
> % Run as is. Nothing needs to be edited.
>
> clear; close all;
>
> %% Check A. Resonance, where the exact answer is known
> % y'' + y = cos(t), y(0) = 0, y'(0) = 0 has exact solution y = (t/2) sin(t).
> m = 1; b = 0; k = 1; f = @(t) cos(t);
> rhs = @(t, Y) [ Y(2); ( f(t) - b*Y(2) - k*Y(1) ) / m ];
> exact = @(t) t/2 .* sin(t);
> tspan = [0 50];
>
> tols = [1e-3, 1e-6, 1e-9]; % 1e-3 is MATLAB's default RelTol
> labels = cell(1, numel(tols));
>
> figure('Position', [100 100 1000 400]);
> subplot(1, 2, 1);
> fprintf('Check A: resonance, compared with y = (t/2) sin t\n');
> fprintf(' RelTol output times max |error|\n');
> for j = 1:numel(tols)
> opts = odeset('RelTol', tols(j), 'AbsTol', tols(j)*1e-3);
> [t, Y] = ode45(rhs, tspan, [0; 0], opts);
> err = abs(Y(:,1) - exact(t));
> fprintf(' %7.0e %10d %10.2e\n', tols(j), numel(t), max(err));
> semilogy(t, err + eps, 'LineWidth', 1.2); hold on;
> labels{j} = sprintf('RelTol = %.0e', tols(j));
> end
> hold off; grid on;
> xlabel('t'); ylabel('|ode45 - exact|');
> ylim([1e-12 1]);
> title('Check A: error against the exact solution');
> legend(labels, 'Location', 'southeast');
>
> %% Check B. Same problem, two tolerances, long run
> % Undamped, unforced: y'' + y = 0, y(0) = 1, y'(0) = 0.
> % The exact orbit is the unit circle, forever.
> m = 1; b = 0; k = 1; f = @(t) 0;
> rhs = @(t, Y) [ Y(2); ( f(t) - b*Y(2) - k*Y(1) ) / m ];
> tspan = [0 1000];
>
> [tL, YL] = ode45(rhs, tspan, [1; 0]); % defaults
> [tT, YT] = ode45(rhs, tspan, [1; 0], odeset('RelTol',1e-8,'AbsTol',1e-10)); % tight
>
> fprintf('\nCheck B: undamped spring on [0, 1000], distance from origin at the end\n');
> fprintf(' default tolerances: %.4f (should be 1)\n', norm(YL(end,:)));
> fprintf(' tight tolerances: %.8f\n', norm(YT(end,:)));
>
> subplot(1, 2, 2);
> plot(tL, sqrt(YL(:,1).^2 + YL(:,2).^2), 'LineWidth', 1.5); hold on;
> plot(tT, sqrt(YT(:,1).^2 + YT(:,2).^2), 'LineWidth', 1.5);
> hold off; grid on; ylim([0.95 1.01]);
> xlabel('t'); ylabel('distance from origin in the phase plane');
> title('Check B: undamped spring, should stay at 1');
> legend('default tolerances', 'RelTol = 1e-8', 'Location', 'southwest');
> ```
> [!example]- MATLAB - `underTheHood.m` (Euler vs. RK4, and where ode45 steps)
> ```matlab
> %% underTheHood.m
> % Two pictures of what a time stepper is doing. Run as is.
> % Part A. Euler and RK4 with the same fixed step, against the exact answer.
> % Part B. The step sizes ode45 actually chooses on a problem with a sudden kick.
>
> clear; close all;
>
> %% Part A. Euler vs. RK4, same step size
> % Underdamped spring: y'' + 0.5 y' + y = 0, y(0) = 1, y'(0) = 0.
> m = 1; b = 0.5; k = 1; f = @(t) 0;
> rhs = @(t, Y) [ Y(2); ( f(t) - b*Y(2) - k*Y(1) ) / m ];
>
> % Exact solution of the underdamped free response
> g = b/(2*m); % decay rate
> wd = sqrt(k/m - g^2); % damped frequency
> y0 = 1; v0 = 0;
> exact = @(t) exp(-g*t) .* ( y0*cos(wd*t) + (v0 + g*y0)/wd * sin(wd*t) );
>
> h = 0.25; % one fixed step size for both methods
> t = 0:h:20;
> N = numel(t);
>
> % Euler: follow the slope at the start of the step.
> YE = zeros(2, N); YE(:,1) = [y0; v0];
> for n = 1:N-1
> YE(:,n+1) = YE(:,n) + h * rhs(t(n), YE(:,n));
> end
>
> % RK4: sample the slope four times across the step and average.
> YR = zeros(2, N); YR(:,1) = [y0; v0];
> for n = 1:N-1
> k1 = rhs(t(n), YR(:,n));
> k2 = rhs(t(n) + h/2, YR(:,n) + h/2*k1);
> k3 = rhs(t(n) + h/2, YR(:,n) + h/2*k2);
> k4 = rhs(t(n) + h, YR(:,n) + h*k3);
> YR(:,n+1) = YR(:,n) + h/6 * (k1 + 2*k2 + 2*k3 + k4);
> end
>
> fprintf('Part A, step h = %g (%d steps)\n', h, N-1);
> fprintf(' Euler max error: %.2e\n', max(abs(YE(1,:) - exact(t))));
> fprintf(' RK4 max error: %.2e\n', max(abs(YR(1,:) - exact(t))));
>
> tf = linspace(0, 20, 1000);
> figure('Position', [100 100 1000 400]);
> subplot(1, 2, 1);
> plot(tf, exact(tf), 'k-', 'LineWidth', 1.5); hold on;
> plot(t, YE(1,:), 'o-', 'MarkerSize', 4);
> plot(t, YR(1,:), 's', 'MarkerSize', 4);
> hold off; grid on;
> xlabel('t'); ylabel('y(t)');
> title(sprintf('Part A: same step h = %g', h));
> legend('exact', 'Euler', 'RK4', 'Location', 'northeast');
>
> %% Part B. Where ode45 puts its steps
> % Lightly damped spring, released from y = 1 and then pushed at t = 10:
> % f(t) = 5 for 10 <= t <= 10.5, and 0 otherwise.
> m = 1; b = 0.2; k = 1;
> f = @(t) 5 * (t >= 10 & t <= 10.5);
> rhs = @(t, Y) [ Y(2); ( f(t) - b*Y(2) - k*Y(1) ) / m ];
>
> % Refine = 1 makes ode45 report only the steps it actually took.
> opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-9, 'Refine', 1);
> [t, Y] = ode45(rhs, [0 30], [1; 0], opts);
> hstep = diff(t);
>
> fprintf('\nPart B: %d steps, smallest %.2e, largest %.2e\n', ...
> numel(hstep), min(hstep), max(hstep));
>
> subplot(2, 2, 2);
> plot(t, Y(:,1), '.-', 'MarkerSize', 6); grid on;
> ylabel('y(t)');
> title('Part B: a push at t = 10 (each dot is one step)');
> subplot(2, 2, 4);
> semilogy(t(1:end-1), hstep, '.-', 'MarkerSize', 6); grid on;
> xlabel('t'); ylabel('step size h');
> ```
> [!example]- MATLAB - `missedPush.m` (a skipped forcing, and MaxStep)
> ```matlab
> %% missedPush.m
> % A short push applied to a spring at rest. With its default settings,
> % ode45 can take a step so long that it jumps over the push and never sees it.
> % Limiting the step size with MaxStep fixes it. Run as is.
>
> clear; close all;
>
> m = 1; b = 0.2; k = 1;
> f = @(t) 5 * (t >= 10 & t <= 10.5); % push of size 5 for half a second
> rhs = @(t, Y) [ Y(2); ( f(t) - b*Y(2) - k*Y(1) ) / m ];
> tspan = [0 30];
> Y0 = [0; 0]; % at rest, so nothing happens before t = 10
>
> optsA = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);
> optsB = odeset('RelTol', 1e-6, 'AbsTol', 1e-9, 'MaxStep', 0.1);
>
> [tA, YA] = ode45(rhs, tspan, Y0, optsA);
> [tB, YB] = ode45(rhs, tspan, Y0, optsB);
>
> fprintf('No MaxStep: largest |y| = %.4f\n', max(abs(YA(:,1))));
> fprintf('MaxStep = 0.1: largest |y| = %.4f\n', max(abs(YB(:,1))));
>
> figure('Position', [100 100 1000 400]);
> plot(tA, YA(:,1), 'o-', 'LineWidth', 1.5, 'MarkerSize', 4); hold on;
> plot(tB, YB(:,1), '-', 'LineWidth', 1.5);
> hold off; grid on;
> xlabel('t'); ylabel('y(t)');
> title('A push at t = 10, starting from rest');
> legend('default step control (missed the push)', 'MaxStep = 0.1', 'Location', 'northwest');
> ```
# Further reading
The course notes these pieces came from, for anyone who wants the derivations. [[(307Su26) Day 16 - Euler's method for systems]] covers rewriting higher-order equations as systems and Euler's method. [[(307Su26) Day 17 - Runge-Kutta methods, RK2 and RK4]] builds RK2 and RK4. [[(307Su26) Day 18 - Consolidating IVP solvers, and adaptive Runge-Kutta-Fehlberg]] covers adaptive stepping with an embedded pair. [[(307Su26) Supplement - Deriving RK4]] derives the RK4 coefficients in full.