[C++] Using RungeKutta order 4 to solve system of ODEsby libelec Tags: c++, diff eq, numerical method 

#1
Dec1212, 09:38 AM

P: 176

(DISCLOSURE: I have already posted this problem in http://math.stackexchange.com/questi...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^{4})). 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^{4} = 16 times). I have the following ODE that describes an object's acceleration when falling from high altitude and being subjected to air friction: ∂^{2}y/∂t^{2} = 9.8 + β e^{y/α} v^{2} Because so far I have only learned to deal with these secondorder 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 = ∂^{2}y/∂t^{2} (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) RungeKutta 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^{2} and g(v, t) = v Now, for RungeKutta 4, I need to calculate 4 constants that will give me an idea of the function increment at the ith 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.
I have the following suspects:
Thanks 



#2
Dec1212, 09:49 AM

P: 2,473

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. 



#3
Dec1212, 10:01 AM

P: 176

I changed the code to reflect the tips you gave regarding the crosstype operations:




#4
Dec1212, 10:21 AM

P: 2,473

[C++] Using RungeKutta order 4 to solve system of ODEs
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. 



#5
Dec1212, 11:38 AM

Mentor
P: 14,433

Here's your problem:
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. 



#6
Dec1212, 01:49 PM

P: 176

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. 



#7
Dec1212, 02:41 PM

Mentor
P: 14,433

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



#8
Dec1212, 03:14 PM

HW Helper
P: 6,189

Hi libelec!
RK does not necessarily give an increase of precision of O(h^{4}). 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^{4}), try a different differential equation. For instance one where the acceleration depends linearly on the position (like a system with a spring in it). 



#9
Dec1212, 03:26 PM

Mentor
P: 14,433

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




#10
Dec1212, 07:27 PM

Engineering
Sci Advisor
HW Helper
Thanks
P: 6,341

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. 



#11
Dec1212, 09:28 PM

HW Helper
P: 6,925

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



Register to reply 
Related Discussions  
Runge Kutta method to solve second order ODE  Calculus & Beyond Homework  1  
2nd Order RungeKutta: 2nd Order Coupled Differential Equations  Calculus & Beyond Homework  3  
Reducing third order ODE to a system of first order equs + 4th order rungekutta  Differential Equations  1  
Reducing third order ODE to a system of first order equs + 4th order rungekutta  Calculus & Beyond Homework  0  
4:th order RungeKutta for system of PDE's  Programming & Computer Science  0 