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.
Best Answer
Your problem might be that you stop the loop the instant that $t=h*i\ge 120$ starting at
i=1
. Since the values ofh
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 computingi*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 fori=0
toN
withdt=h
fori=0
toN-1
and in the last stepi==N
, setdt=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 inh
get compensated.And then you will find that with
h=0.01
andh=0.005
you already exceeded the tradeoff point where the accumulating floating point errors have a greater weight than the discretization errors. Comparingh=0.4
andh=0.2
results in the expected factor of 16.