8  RK4 in Wolfram Language

This chapter translates the classical RK4 method derived in Chapter 7 into executable Wolfram Language. We implement the method from scratch, verify it against the known solution \(e^x\), compare it side by side with Euler’s method, RK2, and RK3, examine why four function evaluations per step deliver such disproportionate accuracy per unit of work, and close by showing how Wolfram’s built-in NDSolve relates to the method we have built.

8.1 Implementing RK4 from Scratch

The classical RK4 method assigns the four slope estimates weights \(1/6\), \(1/3\), \(1/3\), \(1/6\), sampling at the start, two midpoints, and the end of each step. Chapter 7 derived these weights by forcing the Taylor expansion of the numerical increment to match the exact solution through order \(h^4\). Here we translate those coefficients directly into Wolfram Language.

rk4Step[f_, x_, y_, h_] := Module[{k1, k2, k3, k4},
  k1 = f[x, y];
  k2 = f[x + h/2, y + (h/2)*k1];
  k3 = f[x + h/2, y + (h/2)*k2];
  k4 = f[x + h,   y + h*k3];
  y + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
]

The function takes the slope function f, the current state \((x, y)\), and the step size h, and returns the next \(y\) value. To integrate over an interval we iterate rk4Step, carrying both \(x\) and \(y\) forward together as a pair. NestList returns every intermediate state, giving the complete solution trajectory.

rk4Run[f_, x0_, y0_, xEnd_, h_] :=
  NestList[
    {#[[1]] + h, rk4Step[f, #[[1]], #[[2]], h]} &,
    {x0, y0},
    Round[(xEnd - x0)/h]
  ]

We reload the test problem established in Chapter 2.5. Every numerical chapter reloads these definitions so each notebook is self-contained.

f[x_, y_] := y
exactSolution[x_] := Exp[x]
x0 = 0; y0 = 1; xEnd = 2;

Run rk4Run with step size \(h = 0.5\), producing five points on \([0, 2]\), and display each alongside the exact value.

soln = rk4Run[f, x0, y0, xEnd, 0.5];
TableForm[
  Table[
    {soln[[i, 1]], soln[[i, 2]], exactSolution[soln[[i, 1]]],
     Abs[soln[[i, 2]] - exactSolution[soln[[i, 1]]]]},
    {i, 1, Length[soln]}
  ],
  TableHeadings -> {None, {"x", "RK4", "Exact", "Error"}}
]

8.2 Verifying Correctness Against Known Solutions

A single run at one step size confirms the method produces numbers close to the exact solution, but it does not confirm the order. To verify the order we show that halving \(h\) reduces the error by a factor of \(2^4 = 16\).

hValues = {0.5, 0.25, 0.125, 0.0625, 0.03125};
errors = Table[
  Abs[Last[rk4Run[f, x0, y0, xEnd, h]][[2]] - exactSolution[xEnd]],
  {h, hValues}
];

Display \(h\), the RK4 approximation at \(x = 2\), the exact value \(e^2\), and the absolute error.

TableForm[
  MapThread[
    {#1,
     Last[rk4Run[f, x0, y0, xEnd, #1]][[2]],
     exactSolution[xEnd], #2} &,
    {hValues, errors}
  ],
  TableHeadings -> {None, {"h", "RK4 at x=2", "Exact", "Error"}}
]

The ratio of consecutive errors is the most direct check on order.

Table[errors[[i]]/errors[[i + 1]], {i, 1, Length[errors] - 1}]

{13.0129, 14.4237, 15.1898, 15.5892}

Each ratio should be close to \(16\). Small deviations appear at very small \(h\) where floating-point rounding competes with the truncation error. This confirms that the method derived symbolically in Chapter 7 behaves exactly as the derivation predicted.

8.3 Euler, RK2, RK3, and RK4 Side by Side

Chapters 4 through 7 derived each method separately. Here we collect all four into a single comparison. We define each lower-order step function with the same interface as rk4Step and then use a common runner so that only the step function changes between runs.

eulerStep[f_, x_, y_, h_] := y + h*f[x, y]

The Midpoint Method (RK2 with \(b_2 = 1\) and \(c_2 = 1/2\), derived in Chapter 5.6) samples \(f\) once at the midpoint of the interval.

rk2Step[f_, x_, y_, h_] := Module[{k1, k2},
  k1 = f[x, y];
  k2 = f[x + h/2, y + (h/2)*k1];
  y + h*k2
]

Kutta’s 1901 third-order method (derived in Chapter 6.5) uses \(c = \{0, 1/2, 1\}\) and \(b = \{1/6, 2/3, 1/6\}\).

rk3Step[f_, x_, y_, h_] := Module[{k1, k2, k3},
  k1 = f[x, y];
  k2 = f[x + h/2, y + (h/2)*k1];
  k3 = f[x + h,   y - h*k1 + 2*h*k2];
  y + (h/6)*(k1 + 4*k2 + k3)
]

A single runner accepts any step function and returns the complete trajectory as a list of \(\{x, y\}\) pairs.

methodRun[stepFn_, f_, x0_, y0_, xEnd_, h_] :=
  NestList[
    {#[[1]] + h, stepFn[f, #[[1]], #[[2]], h]} &,
    {x0, y0},
    Round[(xEnd - x0)/h]
  ]

Run all four methods at \(h = 0.2\) on \([0, 2]\).

hCmp = 0.2;
solEuler = methodRun[eulerStep, f, x0, y0, xEnd, hCmp];
solRK2   = methodRun[rk2Step,   f, x0, y0, xEnd, hCmp];
solRK3   = methodRun[rk3Step,   f, x0, y0, xEnd, hCmp];
solRK4   = methodRun[rk4Step,   f, x0, y0, xEnd, hCmp];

Plot all four trajectories against the exact solution.

fourFig = Show[
  ListLinePlot[
    {solEuler, solRK2, solRK3, solRK4},
    PlotStyle -> {Red, Orange, Green, Blue},
    PlotLegends -> {"Euler", "RK2", "RK3", "RK4"},
    PlotLabel -> Style["Four Methods at h = 0.2", 13],
    Frame -> True, FrameLabel -> {"x", "y"}],
  Plot[exactSolution[x], {x, x0, xEnd},
    PlotStyle -> {Black, Dashed}]];
Export["08_RK4_in_Wolfram/four-methods.pdf", fourFig];
Export["08_RK4_in_Wolfram/four-methods.svg", fourFig];

Four methods compared at \(h = 0.2\) against the exact \(e^x\) (black dashed).

8.4 Cost vs. Accuracy: Function Evaluations per Step

Each step of each method calls \(f\) a fixed number of times: Euler once, RK2 twice, RK3 three times, RK4 four times. The question is whether the extra evaluations are worth the cost.

TableForm[
  {{"Euler", 1, "O(h)"},
   {"RK2",   2, "O(h^2)"},
   {"RK3",   3, "O(h^3)"},
   {"RK4",   4, "O(h^4)"}},
  TableHeadings -> {None,
    {"Method", "f calls per step", "Order"}}
]

For a fair comparison we hold the total number of \(f\) evaluations fixed at 24 over \([0, 2]\). Euler takes 24 steps of size \(h = 1/12\); RK4 takes 6 steps of size \(h = 1/3\). Both perform exactly 24 evaluations of \(f\), so the comparison measures accuracy per unit of computational work.

totalEvals = 24;
hFor[order_] := (xEnd - x0) * order / totalEvals;

errAt2[stepFn_, h_] :=
  Abs[Last[methodRun[stepFn, f, x0, y0, xEnd, h]][[2]] -
       exactSolution[xEnd]];

TableForm[
  {{"Euler", N[hFor[1]], totalEvals/1, errAt2[eulerStep, hFor[1]]},
   {"RK2",   N[hFor[2]], totalEvals/2, errAt2[rk2Step,   hFor[2]]},
   {"RK3",   N[hFor[3]], totalEvals/3, errAt2[rk3Step,   hFor[3]]},
   {"RK4",   N[hFor[4]], totalEvals/4, errAt2[rk4Step,   hFor[4]]}},
  TableHeadings -> {None,
    {"Method", "h", "Steps", "Error (24 total evals)"}}
]

RK4’s error is many orders of magnitude smaller than Euler’s even though both performed exactly 24 evaluations of \(f\). The four slopes are combined in a carefully weighted average that extracts far more accuracy per evaluation than a single slope ever could.

8.5 Error Convergence: Confirming Fourth-Order Behavior Numerically

The ratio test in the previous section showed one snapshot. A more complete picture comes from plotting the error over a range of step sizes on a log-log scale. A method of order \(p\) satisfies

\[E(h) \approx C \cdot h^p\]

so on a log-log plot the error appears as a straight line with slope \(p\). We compute the error at \(x = 2\) for all four methods over a sequence of step sizes, then both plot the curves and fit the slopes.

hList = Table[2.0^(-k), {k, 1, 8}];
errFn[stepFn_] :=
  Table[
    Abs[Last[methodRun[stepFn, f, x0, y0, xEnd, h]][[2]] -
        exactSolution[xEnd]],
    {h, hList}
  ];
eEuler = errFn[eulerStep];
eRK2   = errFn[rk2Step];
eRK3   = errFn[rk3Step];
eRK4   = errFn[rk4Step];

Plot the four error curves on a log-log scale.

loglogFig = ListLogLogPlot[
  {Transpose[{hList, eEuler}],
   Transpose[{hList, eRK2}],
   Transpose[{hList, eRK3}],
   Transpose[{hList, eRK4}]},
  Joined -> True,
  PlotStyle -> {Red, Orange, Green, Blue},
  PlotLegends -> {"Euler", "RK2", "RK3", "RK4"},
  PlotLabel -> Style["Error at x = 2 vs. step size h", 13],
  Frame -> True, FrameLabel -> {"h", "Error"}];
Export["08_RK4_in_Wolfram/loglog.pdf", loglogFig];
Export["08_RK4_in_Wolfram/loglog.svg", loglogFig];

Log-log plot of error vs. step size. Slopes equal the orders 1, 2, 3, 4.

To read the slopes numerically we fit \(\log(\text{Error}) = p \cdot \log(h) + c\) for each method. The slope \(p\) is the order.

slopeOf[errList_] :=
  LinearModelFit[
    Transpose[{Log[hList], Log[errList]}],
    lgh, lgh
  ]["BestFitParameters"][[2]];

{slopeOf[eEuler], slopeOf[eRK2], slopeOf[eRK3], slopeOf[eRK4]}

{0.91696, 1.93237, 2.93085, 3.92812}

The fitted slopes should be approximately \(1, 2, 3, 4\), confirming that each method achieves exactly the order the derivation promised. The log-log plot makes this visible: the RK4 line is far steeper than the Euler line, meaning RK4’s error shrinks much faster as \(h\) decreases.

8.6 How Wolfram’s NDSolve Relates to What We Just Built

NDSolve is Wolfram Language’s built-in ODE solver. For non-stiff problems it uses an adaptive Runge-Kutta method, specifically the Dormand-Prince pair (RK45), which embeds both a 4th-order and a 5th-order estimate in a single pass of 6 stage evaluations. The two estimates are compared to control the step size automatically: the solver takes smaller steps where the solution changes rapidly.

sol = NDSolve[{y'[x] == y[x], y[0] == 1}, y, {x, 0, 2}];
yND[x_] = y[x] /. sol[[1]];

ndFig = Plot[
  {yND[x], exactSolution[x]},
  {x, x0, xEnd},
  PlotStyle -> {Blue, {Black, Dashed}},
  PlotLegends -> {"NDSolve", "Exact"},
  PlotLabel -> Style["NDSolve vs. Exact Solution", 13],
  Frame -> True, FrameLabel -> {"x", "y"}];
Export["08_RK4_in_Wolfram/ndsolve.pdf", ndFig];
Export["08_RK4_in_Wolfram/ndsolve.svg", ndFig];

NDSolve (blue) tracks the exact solution to within machine precision.

Compare NDSolve’s error at \(x = 2\) with our fixed-step RK4 at \(h = 0.2\).

Abs[yND[2] - exactSolution[2]]

          -7
3.20719 10

Abs[Last[rk4Run[f, x0, y0, xEnd, 0.2]][[2]] - exactSolution[2]]

0.000166857

NDSolve’s error is typically many orders of magnitude smaller. It achieves this not by using a fundamentally different method but by using a smarter strategy: the same weighted-average idea, applied at a step size that adjusts automatically to meet a precision target. Understanding the derivation of RK4 is the prerequisite for understanding why adaptive solvers work and where their guarantees come from.