Take the 2-minute tour ×

The Problem

Use the order 4 Runge-Kutta method to solve the differential equation

$ \frac{\partial^2 y}{\partial t^2} = -g + \beta e^{-y/\alpha }*\left | \frac{\partial y}{\partial t} \right |^{2} $

And corroborate that its global error is O(h^4)

The Mathematical model

I turn the problem into a system of order 1 differential equations:

  • $ \frac{\partial y}{\partial t} = v $
  • $ \frac{\partial v}{\partial t} = -g + \beta e^{-y/\alpha }*\left | v \right |^{2} $

Therefore I define the discretization variables u (for position) and v (for speed) as:

  • v = f(v, u, t)
  • u = g(v, t)

And use the following increments for the Runge-Kutta method of order 4:

For u

  • k1v = h f(vn, un, tn)
  • k2v = h f(vn + 0.5 k1v, un + 0.5 k1u, tn + 0.5 h)
  • k3v = h f(vn + 0.5 k2v, un + 0.5 k2u, tn + 0.5 h)
  • k4v = h f(vn + k3v, un + k3u, tn + h)

For v

  • k1u = h f(vn, tn)
  • k2u = h f(vn + 0.5 k1v, tn + 0.5 h)
  • k3u = h f(vn + 0.5 k2v, tn + 0.5 h)
  • k4u = h f(vn + k3v, tn + h)

And use them in the RK4 expression for each of them:

$ u_{n+1} = u_{n} + \frac{1}{6} (k_{1u} + 2 k_{2u} + 2 k_{3u} + k_{4u}) $

$ v_{n+1} = v_{n} + \frac{1}{6} (k_{1v} + 2 k_{2v} + 2 k_{3v} + k_{4v}) $

NOTE: I first solve for v. To calculate the order of the error, I will solve 120 = h i times with h = 0.1, h = 0.05 and use the result given for h = 0.001 as the "real" value, since I don't know the function that solves the ODE. Then I should corroborate that the absolute value of the "real" minus the result I got from h = 0.1 must be 16 times bigger than what I get when I substract the value I got from h = 0.05 from the "real" value.

The program

I'm using C++ to solve this.

#include <iostream>
#include <math.h>
#include <cmath>
#include <sstream>
#include <fstream>
#include <vector>
#include <cstdlib>

long double rungeKutta(long double h)
    long double alpha = 6629;
    long double beta = 0.0047;

    long double pos = 39068;
    long double speed = 0;

    for (int i = 1; h*i < 120; i++)
        long double k1v = h * (-9.8 + beta * exp(-pos/alpha) * pow(speed, 2));
        long double k1y = h * speed;
        long double k2v = h * (-9.8 + beta * exp(-(pos + 0.5*k1y)/alpha) * pow(speed + 0.5*k1v, 2));
        long double k2y = h * (speed + 0.5*k1v);
        long double k3v = h * (-9.8 + beta * exp(-(pos + 0.5*k2y)/alpha) * pow(speed + 0.5*k2v, 2));
        long double k3y = h * (speed + 0.5*k2v);
        long double k4v = h * (-9.8 + beta * exp(-(pos + k3y)/alpha) * pow(speed  + k3v, 2));
        long double k4y = h * (speed + k3v);

        speed = speed + (k1v + 2.0*(k2v + k3v) + k4v)/6;
        pos = pos + (k1y + 2.0*(k2y + k3y) + k4y)/6;

    return pos;

int _tmain(int argc, _TCHAR* argv[])
    long double errorOne = rungeKutta(0.01);
    long double errorTwo = rungeKutta(0.005);
    long double real = rungeKutta(0.0001);

    cout << fabs(real-errorOne) << endl << fabs(real - errorTwo) << endl;
    return 0;

The results

I find that the error is only reduced by HALF and not to the 1/16th of the first result.

What am I doing wrong?? I've run out of ideas.


share|improve this question
I reproduced that problem. Just to make sure you made no programming mistake. Or we both did the same. –  sonystarmap Dec 11 '12 at 20:38
Thanks. Although that is for the worse, because we still don't know what the problem is :(. –  Heathcliff Dec 11 '12 at 20:54
and above the section "NOTE", you write $ 1/6(k1+k2+k3+k4)$ which should be $ 1/6(k1+2 k2+2 k3+k4)$ like it is written in your code –  sonystarmap Dec 11 '12 at 20:57
is there a context for that equation? –  sonystarmap Dec 11 '12 at 21:01
@macydanim: Yes, it's another writing mistake in the post. Still, in the code you can see I factored that expression as 1/6(k1 + 2.0*(k2 + k3) + k4). The context is an object free-falling from high altitude and being subjected to the forces of gravitational pull (the -9.8 part) and the friction with the air (the function of y), all that having factored out the mass of the object. –  Heathcliff Dec 11 '12 at 21:09
show 1 more comment

3 Answers

You have defined:

  • u' = g(v, t)

but you use it as

  • k1u = h f(vn, un)

instead of

  • k1u = h g(vn, tn)
share|improve this answer
Yes, that was a writing mistake. You can see that I don't commit the same mistake in the other constants, nor in the code. I think. –  Heathcliff Dec 11 '12 at 20:27
So the error must be in the code for the k's. Why don't you use functions to begin with so it matches the RK4 method, and then inline the functions (or let the compiler do it for you). –  ja72 Dec 11 '12 at 20:29
I have, but the same behaviour stands. –  Heathcliff Dec 12 '12 at 13:58
@Heathcliff then update the code in the posting to show it. –  ja72 Dec 12 '12 at 16:35
add comment

Your problem might be that you stop the loop the instant that $t=h*i\ge 120$ starting at i=1. Since the values of h used divide 120 (or any integer), this means that the number of steps performed may be one off from the required number due to numerical errors in computing i*h. This may give an error of size $h$ to the desired end time $120$ that is reflected in an error of size $h$ in the solution values.

So to make absolutely sure that the correct time is integrated, define N=floor(120/h), so that $Nh\le 120<N+1$, loop for i=0 to N with dt=h for i=0 to N-1 and in the last step i==N, set dt=120-N*h.

And indeed, if you follow the actual time during integration by introducing t=0 before the loop and updating t+=h inside the loop, you will find that the loop ends at t==120-h. An alternative to using the actual number of steps in the loop is to change the loop condition to (i-0.5)*h<120, so that rounding errors in h get compensated.

And then you will find that with h=0.01 and h=0.005 you already exceeded the tradeoff point where the accumulating floating point errors have a greater weight than the discretization errors. Comparing h=0.4 and h=0.2 results in the expected factor of 16.

share|improve this answer
add comment

This is an old question of mine, but I'll answer since I see it has had some recent activity. The problem arose due to a misinterpretation of the problem and data given. After I re-ran all the simulations that fed the data to my function, it worked properly.

I don't recall if it was exactly order-4 that I got, so there may be some further problems with my function, but since it was close enough I didn't keep looking.

share|improve this answer
add comment

Your Answer


