- #1
libelec
- 176
- 0
(DISCLOSURE: I have already posted this problem in http://math.stackexchange.com/questions/256393/calculate-runge-kutta-order-4s-order-of-error-experimentally, but found no satisfactory answer)
The problem is this. I need to experimentally check that RK4 method has an error of order 4 (O(h4)). To do this, I have calculated the value the method returns in a certain point for the step values h = 0.1 and h = 0.05 and compared those results to the one I got for the step value h = 0.001. I should check, then, that the difference between h0.001 and h0.1 must be 16 times larger than the difference between h0.001 and h0.05 (since I reduced the step by half, the error must be reduced 24 = 16 times).
I have the following ODE that describes an object's acceleration when falling from high altitude and being subjected to air friction:
∂2y/∂t2 = -9.8 + β e-y/α |v|2
Because so far I have only learned to deal with these second-order equations analytically by turning them into a system of ODEs, I propose the following variable change:
- ∂y/∂t = v(t) (derivative of position with respect to time equals speed)
- ∂v/∂t = ∂2y/∂t2 (derivative of speed with respect to time equals acceleration, or the second derivative of position with respect to time).
Therefore I get the system:
- ∂v/∂t = -9.8 + β e-y/α |v|2 = f(v, y, t)
- ∂y/∂t = v = g(v, t)
Runge-Kutta 4 for this system is expressed then as:
vn+1 = vn + 1/6(k1v + 2k2v + 2k3v + k4v)
yn+1 = yn + 1/6(k1y + 2k2y + 2k3y + k4y)
Where f(v, y, t) = -9.8 + β e-y/α |v|2 and g(v, t) = v
Now, for Runge-Kutta 4, I need to calculate 4 constants that will give me an idea of the function increment at the i-th point. This I have to do for each equation, so I calculate:
- For v:
* k1v = h f(v, y, t)
* k2v = h f(v + 0.5 k1v, y + 0.5 k1y, t + 0.5 h)
* k3v = h f(v + 0.5 k2v, y + 0.5 k2y, t + 0.5 h)
* k4v = h f(v + k3v, y + k3y, t + h)
- For y:
* k1y = h g(v, t)
* k2y = h g(v + 0.5 k1v, t + 0.5 h)
* k3y = h g(v + 0.5 k2v, t + 0.5 h)
* k4y = h g(v + k3v, t + h)
In this cases, I find that the increase by 0.5 h or h in t does nothing, since t doesn't appear in neither function.
I have written the following program to solve this. The problem I have is that, no matter the step sizes I take, I always find that the error is reduced by half when I reduce the step size by half.
Any ideas why this happens?
I have the following suspects:
Thanks
The problem is this. I need to experimentally check that RK4 method has an error of order 4 (O(h4)). To do this, I have calculated the value the method returns in a certain point for the step values h = 0.1 and h = 0.05 and compared those results to the one I got for the step value h = 0.001. I should check, then, that the difference between h0.001 and h0.1 must be 16 times larger than the difference between h0.001 and h0.05 (since I reduced the step by half, the error must be reduced 24 = 16 times).
I have the following ODE that describes an object's acceleration when falling from high altitude and being subjected to air friction:
∂2y/∂t2 = -9.8 + β e-y/α |v|2
Because so far I have only learned to deal with these second-order equations analytically by turning them into a system of ODEs, I propose the following variable change:
- ∂y/∂t = v(t) (derivative of position with respect to time equals speed)
- ∂v/∂t = ∂2y/∂t2 (derivative of speed with respect to time equals acceleration, or the second derivative of position with respect to time).
Therefore I get the system:
- ∂v/∂t = -9.8 + β e-y/α |v|2 = f(v, y, t)
- ∂y/∂t = v = g(v, t)
Runge-Kutta 4 for this system is expressed then as:
vn+1 = vn + 1/6(k1v + 2k2v + 2k3v + k4v)
yn+1 = yn + 1/6(k1y + 2k2y + 2k3y + k4y)
Where f(v, y, t) = -9.8 + β e-y/α |v|2 and g(v, t) = v
Now, for Runge-Kutta 4, I need to calculate 4 constants that will give me an idea of the function increment at the i-th point. This I have to do for each equation, so I calculate:
- For v:
* k1v = h f(v, y, t)
* k2v = h f(v + 0.5 k1v, y + 0.5 k1y, t + 0.5 h)
* k3v = h f(v + 0.5 k2v, y + 0.5 k2y, t + 0.5 h)
* k4v = h f(v + k3v, y + k3y, t + h)
- For y:
* k1y = h g(v, t)
* k2y = h g(v + 0.5 k1v, t + 0.5 h)
* k3y = h g(v + 0.5 k2v, t + 0.5 h)
* k4y = h g(v + k3v, t + h)
In this cases, I find that the increase by 0.5 h or h in t does nothing, since t doesn't appear in neither function.
I have written the following program to solve this. The problem I have is that, no matter the step sizes I take, I always find that the error is reduced by half when I reduce the step size by half.
Code:
#include <math.h>
#include <cmath>
#include <sstream>
#include <fstream>
#include <vector>
#include <cstdlib>
using namespace std;
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 main()
{
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;
}
Any ideas why this happens?
I have the following suspects:
- I have built the system wrong, and that's why no matter what algorithm I use I always get linear reduction
- There's something wrong with the way I calculate the constants
Thanks