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:

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