The next few blogs are applying the Parker-Sochacki method to a pendulum. I am starting with the Simple Pendulum and applied the simulation in Houdini. I will show that a higher order can handle larger time steps.
The simple pendulum is a nonlinear Ordinary Differential Equation but has no damping or driving forces acting on it. I will cover those next.
Here is the equation for the simple pendulum:
This is a Second Order Ordinary Differential Equation where w_0 is a constant. For this example, w_0 is set to 1. To use Parker-Sochacki, this is first broken down into two first order ODEs.
Theta represents angular position in radians and omega represents angular velocity in radians per second. To use Parker-Sochacki, we make sin(theta(t)) an auxiliary variable. Doing this will require a derivative so we have another auxiliary variable for cos(theta(t)).
This requires the chain rule. Let a=sin(theta(t)). Taking the derivative of a leads to
a'(t) = cos(theta(t))*theta'(t)
We already established omega(t) as theta'(t). Here is our list of Ordinary Differential Equations:
Note that I switched a dot and b dot between equal signs. Having this setup makes it possible to turn theta and omega into a Maclaurin Series. Here is the recurrence relationship:
a and b require cauchy products. w does not because w_0^2 is a constant. You can set the initial conditions for theta and omega but note that it is in radians. a(0) = sin(theta(0)) and b(0) = cos(theta(0)).
The final form is a Maclaurin series. That Maclaurin series becomes the initial condition for the next iteration.
This render is my first use of HDRI images and I would like to thank Andrew Maynard for his help on that.
Everything was done in Houdini 15. This first simulation is at the default time step of 1/24 set to order 10 for Parker-Sochacki and took about an hour to render.
The second render is Parker-Sochacki at order 10 doubling the time step. The time step is now every other frame or 2/24:
The pendulum moves twice as fast.
Here is the point of this post: Higher order can handle larger time steps. I multiplied the time step by 10 so that now it is 10/24. Look at a simulation of the pendulum at order 2 and a time step of 10/24:
It quickly crashes. Look at a simulation of the pendulum at order 3 and a time step of 10/24:
It lasts longer than order 2 but it still crashes at taking such large time steps.
Now look at a simulation of the pendulum at order 10 and a time step of 10/24:
It is stable for the entire 240 frames of the simulation. That is the power of using Parker-Sochacki at a higher order.