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

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.

-
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

The problem lies in your conversion to system of first order equations. $y$ is a variable in the RHS but it is not considered in the LHS.

The standard form of Runge-Kutta Mehtod is $\frac{d\mathbf{z}}{dt} = f(\mathbf{z},t)$ for some vector $\mathbf{z}$. In your case, $$\mathbf{z} = (\frac{\partial y}{\partial t},\frac{\partial v}{\partial t})$$

So, $f,g$ must only depend on the elements of $\mathbf{z}$. But the term $\beta e^{-\frac{y}{\alpha}}$ depends on $y$.

UPDATE: Define $\mathbf{z} = (z_1,z_2) = (y,\frac{\partial y}{\partial t})$

Then you have, $$\frac{\partial z_1}{\partial t} = z_2$$ $$\frac{\partial z_2}{\partial t} = -g + \beta e^{-\frac{z_1}{\alpha}} +|z_2|^2$$ There you go!!

-
$z=(y,v)$ . With your notation you would have 2nd order derivatives on the LHS –  sonystarmap Dec 11 '12 at 19:02
But how can I get rid of that y? The only way to get rid of the second derivative is by replacing the derivated variable y by its first derivative, therefore getting the first derivative of a variable. And there's no way of getting rid of y in the RHS of the equation that way. –  Heathcliff Dec 11 '12 at 19:04
REGARDING UPDATE: That's what I am doing right now. In my case z1 = y and z2 = v. –  Heathcliff Dec 11 '12 at 19:24
Yupp. Coded it from scratch in matlab. It comes out to be order 1!! –  dexter04 Dec 11 '12 at 20:05
Can you compare it with results from ode45 ? –  sonystarmap Dec 11 '12 at 20:39

You have defined:

• u' = g(v, t)

but you use it as

• k1u = h f(vn, un)

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