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).

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.

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.

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.)

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).

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.

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).