[C++] Using Runge-Kutta order 4 to solve system of ODEs

[libelec]

(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(h⁴)). 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 h_0.001 and h_0.1 must be 16 times larger than the difference between h_0.001 and h_0.05 (since I reduced the step by half, the error must be reduced 2⁴ = 16 times).

I have the following ODE that describes an object's acceleration when falling from high altitude and being subjected to air friction:

∂²y/∂t² = -9.8 + β e^-y/α |v|²

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 = ∂²y/∂t² (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|² = f(v, y, t)
- ∂y/∂t = v = g(v, t)

Runge-Kutta 4 for this system is expressed then as:

v_n+1 = v_n + 1/6(k_1v + 2k_2v + 2k_3v + k_4v)
y_n+1 = y_n + 1/6(k_1y + 2k_2y + 2k_3y + k_4y)

Where f(v, y, t) = -9.8 + β e^-y/α |v|² 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:

* k_1v = h f(v, y, t)
* k_2v = h f(v + 0.5 k_1v, y + 0.5 k_1y, t + 0.5 h)
* k_3v = h f(v + 0.5 k_2v, y + 0.5 k_2y, t + 0.5 h)
* k_4v = h f(v + k_3v, y + k_3y, t + h)

- For y:

* k_1y = h g(v, t)
* k_2y = h g(v + 0.5 k_1v, t + 0.5 h)
* k_3y = h g(v + 0.5 k_2v, t + 0.5 h)
* k_4y = h g(v + k_3v, 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.

#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

 

[jedishrfu]

You have a lot of doubles listed where you assign them with integer values.

You have some calculations where you're dividing double factors with integer values which may be the actual source of the problem.

Lastly, your for loop uses a terminating comparison that compares a long double to an integer. Its always better to compare apples to apples that way you know its not the compiler.

Also I would put a print statement inside your loop to show intermediate results and to see if the for loop is actually executed that may help you determine where if any your problem lies.

 

[libelec]

I changed the code to reflect the tips you gave regarding the cross-type operations:

long double rungeKutta(long double h)
{
	long double alpha = 6629.0;
	long double beta = 0.0047;

	long double pos = 39068.0;
	long double speed = 0;

	for (double i = 1; h*i < 120.0; 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.0;
        pos = pos + (k1y + 2.0*(k2y + k3y) + k4y)/6.0;

	cout << "i = " << i << ", speed = " << speed << "pos = " << pos << endl;
	}

	return pos;
}

But I get the same results. I added the print statement to check if there is some mistake in the for, but found nothing.

 

[jedishrfu]

Here's a discussion on long double usage that may be helpful to you.

http://www.cplusplus.com/forum/beginner/34088/

Basically, I think you need to use 'L' at the end of your long double constants.

What did the for loop print tell you?

Did the loop cycle as many times as you expected?

Did you see 'k', speed and pos values changing as expected?

I can't comment on your algorithm so I can't help there but checking intermediate results will give you a sense that it is working.

 

[D H]

Here's your problem:

cout << "i = " << i << ", speed = " << speed << "pos = " << pos << endl;

You're using the default precision here. You aren't going to see any change until you get to a step size of 3 seconds or so, and with that, you'll only see a tiny difference in the velocity.

You either need to print with greater precision, or even better, save the results from your smallest step size, compute the difference between the final velocity for step size #i versus step size #0, and print those differences along with the values.

 

[libelec]

Here's your problem:

cout << "i = " << i << ", speed = " << speed << "pos = " << pos << endl;

You're using the default precision here. You aren't going to see any change until you get to a step size of 3 seconds or so, and with that, you'll only see a tiny difference in the velocity.

You either need to print with greater precision, or even better, save the results from your smallest step size, compute the difference between the final velocity for step size #i versus step size #0, and print those differences along with the values.

I don't see what difference it makes. The relation is one is half of the other wheather I print 6 or 10 digits. I mean, if I give it larger precision, it may go from printing 0.06666 and 0.03333 to 0.066661234 and 0.033334564.

EDIT: I think I see what you may be missing: my second post only has the one function rungeKutta(); the whole program is in the first post.

 

[D H]

I see now that you are printing the differences in the final position.

Did you do what jedishrfu suggested in post #4:
- Did you fix your use of double constants such as 0.001?
- Did you verify that the loop cycled as many times as you expected?

There is another problem with your scheme. Any numerical integration technique is subject to two main sources of error: Errors inherent in the technique itself, and errors that result from finite precision arithmetic (e.g., long double). Errors inherent in the technique dominate for large step sizes while errors that result from using finite precision arithmetic dominate for small step sizes. Decreasing the step size *increases* the error when the step size is so small that arithmetic is the dominant source of error.

How do you know that those step sizes of 0.0001, 0.05, and 0.01 seconds are all in the range where technique errors dominate? (Hint: You don't. You might want to rethink your use of step sizes.)

 

[I like Serena]

Hi libelec!

RK does not necessarily give an increase of precision of O(h⁴).
It depends on the differential equation.

To clarify, suppose your beta would be zero (no friction), then there would be no effect of step size at all, since each approximation would be perfect.
That is, RK would effectively be O(1).

To find O(h⁴), try a different differential equation.
For instance one where the acceleration depends linearly on the position (like a system with a spring in it).

 

[D H]

This problem does yield O(h^4) behavior when h is the range where errors inherent in RK4 dominate over numerical errors. No problem will exhibit O(h^4) behavior for all time steps if you use finite precision arithmetic (e.g., IEEE floating point).

 

[AlephZero]

How do you know that those step sizes of 0.0001, 0.05, and 0.01 seconds are all in the range where technique errors dominate? (Hint: You don't. You might want to rethink your use of step sizes.)

An empirical way to do that is keep doubling your step size till the solution goes unstable, and then back off from there. The masochistic way is to do a theoretical stability analyiss for your differential equation!

Based on experience, I would expect the long double issues are a red herring. You shouldn't need to use long double arithmetic to show this effect. Doing everything in double should be good enough.

 

[rcgldr]

I'm not sure what compiler you're using, but note that Microsoft's C / C++ compilers ignore long double and just treat it as double (both are 64 bits). The old 16 bit compilers, Visual C / C++ 2.X or older, were the last ones to support long doubles (80 bit floatling point).