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 #include #include #include #include #include #include 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.