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

Question

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?

Thanks in advance.

Answer 1

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.

Answer 2

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.