2 Notes on Applying Parker-Sochacki to Ordinary Differential Equations

Greetings!

The purpose of this blog post is to discuss what you are supposed to do if you add a constant or multiply by a constant in an Ordinary Differential Equation. I have not seen this explicitly discussed in anything related to Parker-Sochacki and, while it is most likely stated somewhere, have not found it discussed when power series methods are talked about. This is brought up in case this comes up in an equation you need to apply Parker-Sochacki to. I dealt with this issue in 2 simulations I will later go over.

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

Case 1: Multiplying by a constant.

If you have an Ordinary Differential Equation with a constant in front and you are using Parker-Sochacki to get a Maclaurin Series, all you need to do is multiply by that constant at the end.

Example: y’ = 2y

After I discuss the 2nd case, I will show an example

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

Case 2: Adding by a constant

If you have an Ordinary Differential Equation where you add by a constant, during the recurrence relationship, you only add by the constant during the first iteration. I’m not talking about the initial condition. I am talking about the very first iteration. After that, ignore the constant. This might be confusion so I am showing this by example. When I say add by a consant:

Example: y’ = y + 1

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

The example I will demonstrate shows both:

y’ = 2y + 1

Attached is the code:

% Adam Wermus
% May 17, 2015
% This program solves the ODE y’ = 2y+1 with the analytical solution
% y = 1/2(3e^2t – 1) ; 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 = (1/2)*(3*exp(2*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
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
if j == 2
z(2) = 2*z(1) + 1;
end
z(j+1) = 2*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
subplot(2,1,1);
plot(t,y,’r’);
xlabel(‘t’)
ylabel(‘y’)
legend(‘Analytical Solution’);

subplot(2,1,2);
plot(t,zsum,’r’);
xlabel(‘t’)
ylabel(‘y’)
legend(‘Approximate Solution’);

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

—————————————————————————————————————————————————-

Here is an image of the file:

y' = 2y + 1

y’ = 2y + 1

This is the chunk in the code I want to note regarding multiplyng and adding. This part computes the Maclaurin Coefficients:

for i = 2:frames
% Compute Maclaurin Coefficients
for j = 1:order-1
if j == 2
z(2) = 2*z(1) + 1;
end
z(j+1) = 2*z(j)/(j);
end

Remember for Matlab, the index 1 is the initial condition. 0 is the index for the initial condition if this was in another language like C++ or Python. Also, if I was in C++ or Python, I would divide by j+1 (I have the recurrence relationship equation mentioned in the previous blog and don’t want you to be confused because the code looks different),

Notice that I am multiplying by 2:
z(j+1) = 2*z(j)/(j);
That is all you need to do to multiply by a constant. Do what you did before but multiply by a constant.

For adding, I only added by 1 when I was computing the first coefficient(for matlab that is index 2). Everything else, I ignored adding by 1.

if j == 2
z(2) = 2*z(1) + 1;
end

Also, this is just a note on how to handle adding constants. The if statement was not necessary. The for loop could be shifted.

Thank you for going through this. You are awesome.

In the next blog, we will cover circular motion and I will show WHY I want to pioneer Parker-Sochacki into Computer Graphics.

Take Care

-Adam Wermus

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s