Greetings!

This blog post is designed to show how use the Parker-Sochacki Method on 3 Ordinary Differential Equations step by step. Code is attached in Matlab. Images from the code are attached. A rendered Houdini video is also attached for example 2.

The first step I had with Parker-Sochacki was-frankly- to see if I could do it. I tested out the method on Ordinary Differential Equations that had a known exact solution so I knew for sure if the Maclaurin series was getting the same result as the analytical, or exact, solution. This also gave me a lot of practice on how to use the Parker-Sochacki method. It is not a standard formula. Depending on what the ODEs are, they have to be played with to get in a certain form. I like to think of it as driving. It takes some practice. You don’t always brake 500 feet before a stop sign. There are guidelines such as how far away from a car you would like to be or having some judgement on when to brake, accelerate, and by how much.

I tested out Parker-Sochacki on multiple ODEs. The objective was to see if the Maclaurin Series can get the same amount of accuracy as the Exact Solution. On the ODEs I tested…yes it can. This was written in Matlab where I could get a clear plot and data. The three ODEs are worth noting because they cover most of what is encountered with Parker-Sochacki.

1: y’ = y

2: y’ = y^2

3: y’ = sin(t)

The final form of Parker-Sochacki, for these three ODEs are:

where i is the number of coefficients. Generally speaking, the higher i is, the more accurate the solution. t is the time step, which is a constant the user can choose. The beauty with Parker-Sochacki is this: Runge-Kutta 4th Order accuracy is somewhere between fourth and fifth order. So my whole vision in animation simulations is to test the benefit of having higher orders. Is there a benefit to having a Parker-Sochacki solution to cloth when it is set to order 10?

The issue Parker-Sochacki covered is how to find the next coefficients. Here they are for these three ODEs.

1: y’ = y

To compute the next coefficient use the following pattern:

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

You will always have initial conditions. this means that you have an already define y_0.

For this example, I used y_0 = 1

Look at how the next few coefficients are computed:

y_(0+1) = y_0 / (0+1) = 1 / 1 = 1

y_1 = 1

y_(1+1) = y_1 / (1+1) = 1 / 2

y_2 = 1 / 2

This pattern continues. These coefficients are then put into the y(t) formula above. This whole process is repeated for as many times as the user wishes. Attached is the Matlab code that solved this ODE where we have the exact answer and the Parker-Sochacki implementation:

——————————————————————————————————————–

% Adam Wermus

% May 17, 2015

% This program solves the ODE y’ = y with the analytical solution

% y = e^t ; y(0) = 1

% and the Parker Sochacki version

% y(i+1) = y(i) / (i+1)

% Clear everything

clc;

close all;

clear all;

% Analytical solution

t = 0:.1:10; % must be same length as frames

y = exp(t);

% Parker-Sochacki Setup

tau = t(2)-t(1); % time step

order = 100; % order of Maclaurin Series

frames = length(t); % total number of frames

zsum = zeros(1,frames);

% Parker-Sochacki Implementation

z(1) = y(t==0); % this gives the value of y where t=0 (only valid if t=0 just once)

% Power Array

%power_array = zeros(1, order);

% this and the power_array are equivalent

% for i = 1:order

% power_array(i) = tau^(i-1);

% end

power_array = tau.^(0:order-1);

for i = 2:frames

% Compute Maclaurin Coefficients

for j = 1:order-1

z(j+1) = z(j)/(j);

end

% stores the values of the maclaurin series

ysum = 0;

% Compute Power series

for j = 1:order

ysum = ysum + (z(j)*tau^(j-1));

end

zsum(i) = ysum;

% reset initial conditions

z(1) = ysum;

end

% Plot

title(‘dy/dt = y’)

subplot(2,1,1);

%semilogy(t,y,’r’);

plot(t,y,’r’);

title(‘dy/dt = y’);

xlabel(‘t’)

ylabel(‘y’)

legend(‘Analytical Solution’);

subplot(2,1,2);

%semilogy(t,zsum,’r’);

plot(t,zsum,’b’);

title(‘dy/dt = y’);

xlabel(‘t’)

ylabel(‘y’)

legend(‘Approximate Solution’);

% plot(1:frames,zsum); PS method

——————————————————————————————————————–

Here is a plot. This is what happens when this program is run:

NOTE: Matlab’s syntax is different than the standard formula I showed. Matlab’s array index starts at 1. This means that when I said y_0 in my example, if this was in Matlab I need to say y_1. In the code, you can see how the pattern is a little different.

2: y’ = y^2

This ODE is more complicated because it introduces the use of a cauchy product. If you are not familiar with multiplying two series together, note that there is a special formula called a cauchy product to do this:

Attached is Matlab code that solves this ODE. We chose the initial condition, y(0) = -1, where the analytical solution for this is y(t) = 1/(t+1)

——————————————————————————————————————–

% Adam Wermus

% May 22, 2015

% This program solves the ODE y’ = y^2 with the analytical solution

% y = 1/(t+1) ; y(0) = -1

% and the Parker Sochacki version

% y(i+1) = summation from k = 0 to i [y(i)y(k-i)] / (i+1)

% Clear everything

clc;

close all;

clear all;

% Analytical solution

t = 0:.5:10; % must be same length as frames

y = -1./(t+1);

% Parker-Sochacki Setup

tau = t(2)-t(1); % time step

order = 100; % order of Maclaurin Series

frames = length(t); % total number of frames

zsum = zeros(1,frames);

% Parker-Sochacki Implementation

zsum(1) = -1;

z(1) = -1;

% Power Array for time

power_array = tau.^(0:order-1);

for j = 2:frames

% Cauchy product setup

cauchy_product = 0;

% Compute Maclaurin Coefficients

for i = 1:order

cauchy_product = 0;

for k = 1:i

cauchy_product = cauchy_product + (z(k)*z(i+1-k))/i;

end

z(i+1) = cauchy_product;

end

% stores the values of the maclaurin series

ysum = 0;

% Compute Power series

for i = 1:order

ysum = ysum + (z(i)*tau^(i-1));

end

zsum(j) = ysum;

% reset initial conditions

z(1) = ysum;

end

% Plot Exact Solution

subplot(2,1,1);

plot(t,y,’r’);

xlabel(‘t’)

ylabel(‘y’)

legend(‘Analytical Solution’);

% Plot Parker-Sochacki Solution

subplot(2,1,2);

plot(t,zsum,’r’);

xlabel(‘t’)

ylabel(‘y’)

legend(‘Approximate Solution’);

% plot(1:frames,zsum); PS method

——————————————————————————————————————–

Here is a plot when this code is run:

This particular problem was also implemented into Houdini. A bunch of points were placed representing initial conditions and then the Parker-Sochacki solver was made. You will see that some points quickly go down, quickly go up, settle at the asymptote set by this equation, and some are just above the asymptote and coast for a bit before spiking it upwards. This happens quickly so you may need to see this a few times:

3: y’ = sin(t)

This is where Parker-Sochacki is truly different. It introduces the use of auxiliary variables. Substitutions are made to the ODEs to get it in the necessary form to compute all of the coefficients. It takes some practice. This will be easier to show.

Let a = sin(t)

so

y’ = a

Now take the derivative of a.

a’ = cos(t)

Let b = cos(t)

b’ = -sin(t)

But OOH! -sin(t) is just -a. So

b’ = -sin(t) = -a

We have the following equations that depend on each other so we can soon setup a recurrence relationship:

y’ = a

a’ = b

b’ = -a

Give it some initial conditions and we can use the following recurrence relationship:

For this, y_0 initial condition can be whatever you choose. a_0 and b_0 are dependent on the equation. For the code you will see in a moment, y_0 is set to -1. a_0 = sin(0)=0, and b_0 = cos(0) = 1. Here is the code:

——————————————————————————————————————–

% Adam Wermus

% June 5, 2015

% This program solves the ODE y'(t) = sin(t)

% The analytical solution is y(t) = -cos(t)

%clear everything

close all;

clear all;

clc;

% Actual Solution

t = 0:pi/20:2*pi;

y = -cos(t);

% Plot Actual Solution

subplot(2,1,1);

plot(t,y,’g’)

title(‘Actual Solution’);

xlabel(‘t’);

ylabel(‘cos(t)’);

% Declare variables for Parker-Sochacki Implementation

order = 100; % order for power series

tau = t(2) – t(1); % time step

frames = length(t); % number of iterations

% Declare a and b

a = zeros(1, order);

b = zeros(1, order);

z = zeros(1, order);

% Initial Conditions

z(1) = -1;

a(1) = 0; % a = sin(t)

b(1) = 1; % b = cos(t)

% Declare a_sum and b_sum. These are the power series sums for a and b

z_sum = zeros(1,frames);

a_sum = zeros(1,frames);

b_sum = zeros(1,frames);

% Give a_sum and b_sum initial values

z_sum(1) = -1;

a_sum(1) = 0;

b_sum(1) = 1;

% Number of frames to iterate through

for iterate = 2:frames

% Calculate Coefficients

for i = 1:order-1

z(i+1) = a(i)/(i);

a(i+1) = b(i)/(i);

b(i+1) = -a(i)/(i);

end

% Create temporary sums for a and b

z_temp_sum = 0;

a_temp_sum = 0;

b_temp_sum = 0;

% Create Maclaurin Series

for i = 1:order

z_temp_sum = z_temp_sum + z(i)*tau^(i-1);

a_temp_sum = a_temp_sum + a(i)*tau^(i-1);

b_temp_sum = b_temp_sum + b(i)*tau^(i-1);

end

% Set a_sum and b_sum to the temp sums

z_sum(iterate) = z_temp_sum;

a_sum(iterate) = a_temp_sum;

b_sum(iterate) = b_temp_sum;

% Reset Initial Conditions

z(1) = z_temp_sum;

a(1) = a_temp_sum;

b(1) = b_temp_sum;

end

% Plot Parker-Sochacki Implementation

subplot(2,1,2);

plot(t,z_sum,’g’);

title(‘Parker-Sochacki Solution’);

xlabel(‘t’);

ylabel(‘Approximate Solution’);

——————————————————————————————————————–

Here is the attached image:

This is confusing at the first read and this is not the only way to do this problem. Parker-Sochacki takes practice. The idea is to get a system of equations that depend on each other in order to get it into a form where you can compute all of the coefficients that you want.

Take Care!

Adam Wermus