# How to implement Heun's method to solve a 2nd order ODE? [closed]

I am trying to simulate a soft body object by having a collection of points connected by springs. I also want to have the object bounce on a plane. I have been able to implement Euler's method to make this simulation work, but it is not satisfactory because when I increase the spring constants it becomes unstable. Ideally I would like to implement the Runge-Kutta 4th order method for this simulation, but for now I just want to implement Heun's method, which is the Runge-Kutta 2nd order method.

From what I understand, this is what you have to do:

1. Compute the acceleration of every vertex.
2. Using Euler's method, compute the new velocities and positions of each vertex after one time interval, and store these new vertices in a seperate data structure from the original vertices.
3. Compute the acceleration of every vertex in their new positions.
4. Compute the average of the 2 accelerations you computed.
5. Using this average acceleration, use Euler's method again to compute the new velocites and positions of each vertex after one time interval, starting from the original vertex velocities and positions.

This almost works, but when the object bounces off the plane, it only bounces half as high each time, essentially losing energy. If I just use Euler's method, it does not lose energy.

I factor the collision and gravity forces when I compute the acceleration. The way I am calculating the collision force from the ground is if a vertex is below the ground, then I pretend there is a spring with an equilibrium distance of 0 between the ground and the vertex. Is there something in my algorithm that is incorrect?

-
There are better ODE-solving algorithms than RK when stiffness is a concern. My hands-down favorite, if the system is linear, as it sounds like yours is, is matrix-exponent. – Mike Dunlavey Dec 4 '11 at 14:50
"This almost works, but when the object bounces off the plane, it only bounces half as high each time, essentially losing energy." One possible cause of that is transfer of kinetic energy from overall body motion to internal motion. – mmc Dec 4 '11 at 19:02
@mmc: if this were the answer, it would be the same in Euler. – Ron Maimon Dec 7 '11 at 7:07
@RonMaimon Maybe he is having gross energy non-conservation when using Euler, but you are probably right. – mmc Dec 7 '11 at 11:51
@Thomas Ryabin: Maybe the question would fit better on scicomp.SE? Don't forget to provide crosspost links if you do this (without the help of moderators), i.e. crossposting (as opposed to migration). Alternatively, you can flag the moderators for migration. – Qmechanic Mar 4 '12 at 16:15

## closed as off topic by Qmechanic♦Dec 30 '12 at 7:46

Questions on Physics Stack Exchange are expected to relate to physics within the scope defined in the FAQ. Consider editing the question or leaving comments for improvement if you believe the question can be reworded to fit within the scope. Read more about closed questions here.

Your description of the algorithm is correct, and this is the best quick-and-dirty algorithm for day-to-day simulations. But you have to specify far more detail if you want a certain answer.

One possible reason is that the force near the ground is too sharply varying. You can replace the ground-particle force with just reflecting boundary conditions, so that whenever the y-coordinate of any particle is less than zero, you reverse it's y-velocity (instantly, without any Heun back and forth, just sweep over the particles at the beginning of each time step, and if y<0, $v_y=-v_y$). Then you don't need to have any force from the ground at all. But you will get sound losses in the material during the bounce.

If you want a good diagnostic, add up the total energy (potential + kinetic) and plot it with time. It is difficult to say more, because you can't debug code without seeing it.

-
 The way I am calculating the force from the ground is if a vertex is below the ground, then I pretend there is a spring with an equilibrium distance of 0 between the ground and the vertex. I'll try putting the collision calculations outside of the Heun algorithm. – Thomas Dec 4 '11 at 17:26 The half spring should work if the spring constant is not too big. – Ron Maimon Dec 4 '11 at 17:43

If you have the accelerations as a=f(t,x,v) then starting from t_1, x_1, v_1 and taking a time step h, do the following:

1. K0 = f(t_1, x_1, v_1)
2. K1 = f(t_1+h, x_1+h*v_1, v+h*K0)
3. t_2 = t_1 + h
4. v_2 = v_1 + h*(K0/2+K1/2)
5. x_2 = x_1 + h*(v_1+h*(K0/3+K1/6))

where x_2 and v_2 are the positions and velocities after the step h.

-
 Where did the "K0/3 + K1/6" in step 5 come from? – Thomas Dec 5 '11 at 19:03 If you assume a linear variation of acceleration with respect to time it comes out like that. The standard method has a constant acceleration x_1+h*(v_1+h*K0/2) term only and it is way too unstable. – ja72 Dec 5 '11 at 22:58 @Ja72: It is only "unstable" in the technical sense of stiffness, which happens when you take time scales which are larger than the largest inverse frequency in the problem. This is never the case in discrete particle/spring simulations. He could use the obvious Euler algorithm, or any of the simplest correctors and the answer should be grossly the same. It is pointless to make a better algorithm for this type of thing. – Ron Maimon Jan 4 '12 at 15:03