I wrote an RK4 algorithm and am testing it on $y' = -ty$ which has the solution $y(t) = e^{-t^2 / 2}$

I decided to also graph the error, which I am now trying to decipher. I plotted the solution and error for several time steps. I halved the time step for each test:

$\Delta t = 0.1$ : Delta T = 0.1

$\Delta t = 0.05$ : enter image description here

$\Delta t = 0.025$ : enter image description here

In general, my questions are:

How can the error in my algorithm be analyzed? And, why does that first graph look so funny?

As far as the work I've done:

I am conceptually aware of the following jargon and roughly what it means, however the class I took in numerical methods was run by a genius, not a communicator:

  • "truncation error"
  • "round-off error"
  • RK4 is "fourth order" and so error should drop like $(\Delta t)^4$

My best guess is that error increases with each time step due to round-off. However, I have no idea how to assure myself that this is true.

The only thing I knew to do was to check that error is proportional to $(\Delta t)^4$. It sure seems like the error isn't dropping that fast. I found the maximum error for the first test:

$E_{\Delta t} = c(0.1)^4 = 0.4321$

And found c = 4321, and applied it to

$E_{\Delta t /2} = 0.0930$

but $4321(0.05)^4 = 0.027$, which is roughly a third of the error that I got. I remember my professor mentioning "order of magnitude" a lot. I guess they are within an "order of magnitude", so does that mean everything is good here?

def rk4(dt, t, field, y_n):

    k1 = dt * field(t, y_n)
    k2 = dt * field(t + 0.5 * dt, y_n + 0.5 * k1)
    k3 = dt * field(t + 0.5 * dt, y_n + 0.5 * k2)
    k4 = dt * field(t + 0.5 * dt, y_n + k3)

    return y_n + (k1 + 2 * k2 + 2 * k3 + k4) / 6

if __name__ == '__main__':

    # the ODE y' = -t * y, which has solution y = exp(-t^2 / 2)
    def field(t, vect):
        return np.array([-t * vect])

    # Set the interval over which we want a solution.
    t_0 = -10
    t_n = 10
    dt = .05

    # Determine number of steps in accordance with mesh size
    steps = int((t_n - t_0) / dt)
    time = np.linspace(t_0, t_n, steps, endpoint=False)
    # time = np.arange(t_0, t_n, dt)

    # Initialize solution vectors and error collection
    x = np.zeros(steps)
    error = np.zeros(steps)
    x[0] = 1.928749848e-22
    error[0] = 0

    for i in range(1, steps):
        x[i] = rk.rk4(dt, time[i-1], field, x[i-1])
        error[i] = abs(x[i] - math.pow(math.e, (-time[i] ** 2) / 2)) / math.pow(math.e, (-time[i] ** 2) / 2)
  • 1
    $\begingroup$ Please check your algorithm carefully. From second to third the error curve is only divided by two, not by about 16 as the RK4 method would have you expect. This order-one behavior hints to some error, wrong constant, wrong sign,... $\endgroup$ Commented Oct 16, 2019 at 5:58
  • 1
    $\begingroup$ Another possibility is that you have an indexing error so that the numerical solution and exact solution are evaluated at points with distance $Δt$. Then the resulting error of the shift, $(y(t+Δt)-y(t))=O(Δt)$, would dominate the $O(t(Δt)^4)$ method error. $\endgroup$ Commented Oct 16, 2019 at 6:05
  • $\begingroup$ @LutzL I will do both of your suggestions. I enjoyed learning a little bit from your answer below, as far as being more concise with my code (especially in the realm of graphing). I added mine to the question, so you can see how ugly it is. $\endgroup$ Commented Oct 16, 2019 at 15:41
  • $\begingroup$ In k4 you forgot to remove one of the copied 0.5*. This is sufficient to destroy the order conditions and reduce the order to 1. $\endgroup$ Commented Oct 16, 2019 at 22:40

1 Answer 1


Implementing the RK4 method as

def RK4integrate(f,t,y0):
    y = np.asarray(len(t)*[y0]);
    for i in range(len(t)-1):
        h = t[i+1]-t[i];
    return y

and producing a combined plot of solution graphs and error profiles for the relative error divided by the expected scale $h^4$ by

def p(t): return np.exp(-t**2/2)
def odefunc(t,x): return -t*x 

fig, ax = plt.subplots(2,1,figsize=(12,10))
t0, tmax=-10, 10
for h in [0.1, 0.05, 0.025, 0.01, 0.005 ][::-1]:
    t = np.arange(t0,tmax,h);
    y = RK4integrate(odefunc, t, np.array([p(t[0])]));
    ax[0].plot(t,y[:,0],'-o', ms=1+13*h, label="h=%.3g"%h);
    ax[1].plot(t,(y[:,0]/p(t)-1)/h**4,'-o', ms=1+16*h, label="h=%.3g"%h);
for gr in ax: gr.grid(); gr.legend();

produces the plot

plot of numerical solutions and error profiles

where the convergence of the error profiles shows clearly that the method has order 4 and that the transition from $e^{-50}$ at $t_0=-10$ to the value $1$ at $t=0$ does produce a relatively benign relative error of about $1500\,h^4$ at $t=0$.

Additionally, with smaller step sizes the error profile more and more reflects the symmetry of the problem, meaning that the errors at $t>0$ have the opposite sign but about the same size as the error at $-t$ so that they compensate. This means that the error coefficient at $t=10$ is zero for $h^4$ and what can be seen is the $h^5$ term, accounting for the halving in the scaled relative error at each step size halving.

  h      relative error          scaled rel. error
0.005   5.9285699682831705e-08    94.85711949253073 
0.01    1.8950046616339478e-06    189.50046616339478 
0.025   0.00018492185995810928    473.39996149275964 
0.05    0.005975343139402733      956.0549023044372 
0.1     0.21902043404195348      2190.204340419534 

Away from that point the errors behave as expected for a fourth order method.

Conclusion: Your observed error curve is not reproducible.

  • $\begingroup$ And may I also ask about your relative error calculation? I was taught that relative error is (actual - estimate) / actual. However, I see you did (estimate / actual / stepsize). What is the reason for this? $\endgroup$ Commented Oct 16, 2019 at 15:52
  • $\begingroup$ Having reworked my code to the extent I thought appropriate without totally plagiarizing yours, I still get the exact same error curves. I checked all indices to ensure no 'off by one' errors. In fact, I converted my code from linspace to arange to verify that my indices matched yours. The only difference I noted was that y[i+1,:]=y[i]+(k1+2*k2+2*k3+k4)/6;, and am curious why you used slicing here, as it seems that y[i + 1] would suffice. Indeed, any other way throws an error in my code. $\endgroup$ Commented Oct 16, 2019 at 19:13
  • $\begingroup$ I'm not really sure what it actually does, the idea was to overwrite the existing array instead of replacing one reference by another. In the balance, there would be no real difference in the memory management, one of the arrays in the assignment will need to be destroyed afterwards. $\endgroup$ Commented Oct 16, 2019 at 22:46

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .