Calculate Runge-Kutta order 4's order of error experimentally

Question

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;
    system("pause");
    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.

Thanks.

Answer 1

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.

Answer 2

You have defined:

• u' = g(v, t)

but you use it as

• k1u = h f(vn, un)

instead of

• k1u = h g(vn, tn)

Answer 3

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.