This post is designed to be an introduction to the research I do. This post compares the Parker-Sochacki solution to the motion of a circle with Forward Euler, Runge-Kutta 2nd Order Midpoint Method, Runge-Kutta 3rd Order, and Runge-Kutta 4th Order in terms of time steps. I have code in Matlab, step by step instructions for different numerical integration methods, and videos in Houdini. When I was an intern at Side Effects in the fall of 2015, I created an example file entitled “CircleSolvers” as an educational tool. If you are inside Houdini’s help, copy and paste this link. It will take you straight to it:

http://127.0.0.1:48626/examples/nodes/sop/solver/CircleSolvers

x'(t) = y(t) ; x(0) = 1

y'(t) = -x(t) ; y(0) = 0

To use Parker-Sochacki, first create the Maclaurin Coefficients. We have initial conditions so you can compute as many coefficients as you want.

x_(i+1) = y_i / (i+1)

y_(i+1) = -x_i / (i+1)

Then put it in the following form:

x(t) = x_0*t^0 + x_1*t^1 + x_2*t^2 + x_3*t^3 + … x_n*t^n

y(t) = y_0*t^0 + y_1*t^1 + y_2*t^2 + y_3*t^3 + … y_n*t^n

where t is the time step, a constant you choose. Also, t^0 is 1. I wrote t^0 in case the reader is not familiar with a Taylor Series or a Maclaurin Series and needs to see the pattern.

This first test was implementing it into Matlab. Here is the code and a perfect circle:

Matlab Code:

————————————————————————————————————————————————————

%Adam Wermus

%8/19/14

%This uses the Parker-Sochacki Method as a numerical method to trace a

%circle

%

% x'(t) = y(t) ; x(0) = 1

% y'(t) = -x(t) ; y(0)= 0

% There are two equations in our system

% Setup t and N

N=160;

t=0.01;

% Initial Conditions and number of points

r=1;

s=0;

numPts = 2000;

% answerX and answerY are being set to zero.

answerX=zeros(1,numPts);

answerY=zeros(1,numPts);

% Set up iterative process

for j=1:numPts

% Initial Values for the two IVODEs

x(1)=r;

y(1)=s;

% Iterative Process to find coefficients

for i=1:N

x(i+1)=(y(i) / (i));

y(i+1)=(-x(i) / (i));

end

% Intialize r and s to zero

r=0;

s=0;

% Power Series Solution answerX and answerY

for i=1:N

r=r+x(i)*t^(i-1);

s=s+y(i)*t^(i-1);

end

% This becomes the next term in the summation

answerX(j)=r;

answerY(j)=s;

end

% Plot a circle and our PS method to compare

angle=linspace(0,2*pi);

plot(answerX, answerY, sin(angle),cos(angle), ‘r’);

————————————————————————————————————————————————————

Here is the Plot:

When a numerical integration method is not accurate for this example, rather than trace a circle, it will diverge spiraling inward or outward. Generally speaking, taking a smaller time step will increase the accuracy. However, sometimes that isn’t good enough.

To see what I mean, look at the following video I made of a point moving in a perfect circle using the exact answer. There is another point there moving in a circle moving at a faster rate. This is using the Forward Euler method. Because it is not accurate, by the time the point meets the other point, it is already diverging away from circular motion:

Euler Circle:

Compare that with Parker-Sochacki, at a higher order than 1, which traces the circle perfectly.

Parker-Sochacki Circle:

_________________________________________________________________________________________________

This is a tangent. In case you are learning numerical methods: Here is how I solved the motion for the two Ordinary Differential Equations that solve the motion of a circle. I am copying pages from my Master’s Thesis:

_________________________________________________________________________________________________

Note that “TimeInc” represents the time step. Earlier I called it “t.”

Back to the point of this blog. Inside Houdini, Parker-Sochacki was written at order 15. In terms of time steps, it out performed the other methods. I did a test to see, over 2400 frames(this was meant to be long), when could I break the other methods?

Forward Euler: No matter how small I made the time step, even 0.001, it still diverged from the circle quickly. There is not point discussing Forward Euler any further for this example.

The default time step for Houdini is 1/24 or about 0.04167. This means that there are 24 frames in 1 second. All of the other methods were stable under this time step. How much larger could I take the time step until the particle no longer traced a circle?

Runge-Kutta 2nd Order: I was able to increase the time step to 0.08 before the particle started to spiral outwards

Runge-Kutta 3rd Order: I was able to increase the time step to 0.12 before the particle started to spiral outwards

Runge-Kutta 4th Order: I was able to increase the time step to 0.30 before the particle started to spiral outwards.

Parker-Sochacki 15th Order: I was able to increase the time step to 1.1 before the particle started to spiral outwards.

This is what happens when you take a larger time step but can still maintain an accurate simulation:

A larger time step means that the simulation does not need to be recomputed as many times.

Another test was to see how long it takes to compute the simulation. This is such a fast example that all of the methods were around 0.16-0.17 seconds. This is promising because Parker-Sochacki held its own.

Now this is why all of this is important. I could extend Parker-Sochacki to a higher order and make the time step even larger. It was the best choice if I had to solve for circular motion next to using the analytical solution. Most Ordinary Differential Equations do not have a known analytical solution. Its accuracy is greater than the other methods used because of its higher order. This is just circular motion, however, it’s possible to apply this to any physical simulation with Ordinary Differential Equations.