Recurrence Relation Note

Hello, this is to address a note on recurrence relations and why there have been two cases now where I computed an “i” term rather than an “i+1” term.

Here I will repeat the Projectile with Gravity and Drag and the 2-Body problem recurrence relationship

Projectile with Gravity and Drag Forces:

Baseball ODEs6

2D-2 Body Problem:

PS_2_body_recurrence

For the projectile, p_i and q_i are the terms computed with just “i” and that goes for s2.s3, and a in the  in the 2D-2 Body problem.

 

When setting these equations up, you cannot multiply more than 2 variables at a time. The cauchy product is designed to multiply 2 variables at a time. If there are more(in both cases above there are), another variable is there as a placeholder. In the 2D 2-Body problem, “s” is computed by multiplying “a” and what would be “s^3.” To do this, an “s2” variable is made that is a cauchy product of s and s. Then an “s3” variable is made that is a cauchy product of “s2” and “s” to update “s”  completing the recurrence setup. Also, the variable”a” is multiplying a difference, which can be done. The velocity variables are doing that too.

I’ll demonstrate an example that shows this idea.

y’ = y^3

This has a known analytical solution of 1/sqrt((1-2t)). There is a specified domain so that it doesn’t blow up.

A variable needs to be created that is a cauchy product of two y terms so that the final way to update y is multiplying two variables.

let a = y^2

so

y’ = a*y

Now here is the recurrence relationship:

yprimeequalsycubedr.PNG

We start with a y(0) that is the given initial condition. Then that is used to compute a(0). This cycle repeats and then you use that to get the Maclaurin Series for only y(t). This is a placeholder and so the initial condition for the next frame or calculation is not the Maclaurin Series.

Below is Matlab Code:

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

% Adam Wermus
% November 16, 2015

% This program solves the ODE y’ = y^3 with the analytical solution
% y = 1/sqrt((1-2t)) ; y(0) = 1
% and the Parker Sochacki version
% Clear everything
clc;
close all;
clear all;
% Analytical solution
t = -10:.01:-5; % must be same length as frames
z = 1./sqrt(1-(2.*t));

% Parker-Sochacki Setup
tau = t(2)-t(1); % time step
order = 100; % order of Maclaurin Series
frames = length(t); % total number of frames
y_sum = zeros(1,frames);
a_sum = zeros(1,frames);
a = zeros(1,order);
y = zeros(1,order);

% Parker-Sochacki Implementation
y_sum(1) = 0.2182;
y(1) = y_sum(1);
a_sum(1) = y_sum(1).*y_sum(1);
a(1) = a_sum(1);

% Power Array for time
power_array = tau.^(0:order-1);

for j = 2:frames

% Cauchy product setup
cauchy_product_y = 0;
cauchy_product_a = 0;

% Compute Maclaurin Coefficients
for i = 1:order
cauchy_product_y = 0;
cauchy_product_a = 0;
for k = 1:i
cauchy_product_y = cauchy_product_y + (a(k)*y(i+1-k))/i;
cauchy_product_a = cauchy_product_a + (y(k)*y(i+1-k));
end
y(i+1) = cauchy_product_y;
a(i) = cauchy_product_a;
end

% stores the values of the maclaurin series
y_temp_sum = 0;
a_temp_sum = 0;

% Compute Power series
for i = 1:order
y_temp_sum = y_temp_sum + (y(i)*tau^(i-1));
a_temp_sum = a_temp_sum + (a(i)*tau^(i-1));
end
y_sum(j) = y_temp_sum;
a_sum(j) = a_temp_sum;

% reset initial conditions
y(1) = y_temp_sum;
a(1) = a_temp_sum;
end
% Plot Exact Solution
subplot(2,1,1);
plot(t,z,’r’);
title(‘dy/dt = y^3’)
xlabel(‘t’)
ylabel(‘y’)
legend(‘Analytical Solution’);

% Plot Parker-Sochacki Solution
subplot(2,1,2);
plot(t,y_sum,’r’);
title(‘dy/dt = y^3’)
xlabel(‘t’)
ylabel(‘y’)
legend(‘Approximate Solution’);

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

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

Here is a plot with the analytical solution on top and the Parker-Sochacki solution on the bottom.

yprimeequalsycubed.PNG

 

I hope that clears up why some recurrence relations have only an i term in them.

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