Parker-Sochacki Conservation of Energy on the Simple Pendulum

Back to the Parker-Sochacki Simple-Pendulum: Conservation of Energy

In previous posts, I applied the Parker-Sochacki method to three types of Pendulums: Simple Pendulum, Driven Pendulum, and a Damped Driven Pendulum. I demonstrated how the Parker-Sochacki method uses auxiliary variables to create a Maclaurin series for the angular displacement and angular velocity up to any order. When the order of the polynomial that is the Maclaurin series was raised, larger time steps could be taken without the pendulum “blowing up.” I showed the math and simulations rendered in Houdini. These posts were done towards the end of 2015 if you want to see that section.

Up until now, I have showed posts applying Parker-Sochacki to Initial Value Ordinary Differential Equations, a particle tracing a circle, projectile with gravity and drag forces, a section of a fluid simulation, a section of a wind drag simulation, the two-body/three-body/n-body gravity problem, a simple/driven/damped driven pendulum, and a single/double 1D and 2D mass-spring showing math, code, and rendered animations in Houdini. I have not discussed the topic of energy and thought it would be good to show experimentally if the Parker-Sochacki method conserves energy when applied to the Simple Pendulum. I will show code and results under different orders and time steps.

For this blog, I focus on energy and not how Parker-Sochacki is applied. However, you can see this in the code. If you want to see the math step-by-step, refer to the following post:

(scroll down as it starts on the driven pendulum)

This simulation is a Simple Pendulum. This means that wherever I start the pendulum, in theory, it should swing back and forth without damping or gaining energy. If that happens, that is caused by numerical dissipation. The total energy of the system is calculated by summing the Potential Energy and Kinetic Energy. Parker-Sochacki creates a Maclaurin series for angular displacement of the pendulum bob theta and angular velocity w. To calculate the Potential Energy and Kinetic Energy:

TE = PE + KE

h = length * (1 – cos(theta))

v = length * w

PE = mass * g * h

KE = 0.5 * mass * v * v

Let’s break that down. Length refers to the length of the pendulum rod. h is calculating a “height” for the pendulum. g refers to acceleration due to gravity, which is a constant. We are calculating the angular velocity of the pendulum, which means that if I want to calculate the Kinetic Energy of the pendulum, then I need to convert the angular velocity to linear velocity v., which is multiplying the angular velocity and length of the pendulum.

I will show the code now so you can see what values were chosen. Every “frame,” the Potential Energy PE and Kinetic Energy KE are calculated. Here is the code, followed by a plot of Theta vs. Time


% Adam Wermus
% September 28, 2016
% This program solves the motion of a simple pendulum
% theta”(t) = -w_0^2*sin(theta(t))

% List of First Order Ordinary Differential Equations
% theta'(t) = w(t)
% w'(t) = -w_0^2*a(t)
% a'(t) = b(t)*w(t)
% b'(t) = -a(t)*w(t)
% clear everything
clear all;

% Setup constants
g = 10; % gravity
length = 1; % length of pendulum rod
mass = 1;
constant = g/length;
% Setup Time Step and Order
Order = 20;
multiplier = 1;
Time_Step = multiplier.*(0.041671);
idea = 1:Time_Step:24;
Frames = numel(idea);

% setup array for recurrence
iterate = zeros(1, Order);

% Compute division part of iterations so this is done once
for i = 1 : Order

iterate(i) = 1./i;


% setup array for time
Time = zeros(1, Order);

% Compute Time series
for i = 1: Order

Time(i) = Time_Step^(i-1);


% Setup arrays for theta, w, and auxiliary variables
theta = zeros(1, Order);
w = zeros(1, Order);
a = zeros(1, Order);
b = zeros(1, Order);

% Setup initial conditions
theta(1) = pi/2;
w(1) = 1;
a(1) = sin(theta(1));
b(1) = cos(theta(1));

% Setup Sums
theta_sum = zeros(1, Frames);
w_sum = zeros(1, Frames);
a_sum = zeros(1, Frames);
b_sum = zeros(1, Frames);

% Setup Initial Sums
theta_sum(1) = theta(1);
w_sum(1) = w(1);
a_sum(1) = a(1);
b_sum(1) = b(1);

% Setup Potential, Kinetic, and Total Energies
Potential = zeros(1, Frames);
Kinetic = zeros(1, Frames);
Total_Energy = zeros(1, Frames);

% First Potential Energy
h = (1 – cos(theta(1)));
h = length * h;
Potential(1) = mass*g*h;

% First Kinetic Energy
v = length * w(1);
Kinetic(1) = 0.5*mass*v*v;
% First Total Energy
Total_Energy(1) = Potential(1) + Kinetic(1);

% Go through all frames
for j = 2:Frames

for i = 1:Order-1

% Update Theta Coefficients
theta(i+1) = iterate(i).*w(i);

% Update w Coefficients
w(i+1) = -constant*iterate(i).*a(i);

% Cauchy Product Section for auxiliary variables a and b
cp_a = 0;
cp_b = 0;

for k = 1:i

cp_a = cp_a + (b(k).* w(i-k+1));
cp_b = cp_b + (a(k).* w(i-k+1));


% Update w, a, and b
cp_b = -cp_b;
a(i+1) = iterate(i).*cp_a;
b(i+1) = iterate(i).*cp_b;


% Setup temporary sums
theta_temp_sum = 0;
w_temp_sum = 0;
a_temp_sum = 0;
b_temp_sum = 0;

% Compute Maclaurin Series
for i = 1 : Order-1

theta_temp_sum = theta_temp_sum + (theta(i) * Time(i));
w_temp_sum = w_temp_sum + (w(i) * Time(i));
a_temp_sum = a_temp_sum + (a(i) * Time(i));
b_temp_sum = b_temp_sum + (b(i) * Time(i));


% Store Values into sum
theta_sum(j) = theta_temp_sum;
w_sum(j) = w_temp_sum;
a_sum(j) = a_temp_sum;
b_sum(j) = b_temp_sum;

% Reset Initial Conditions
theta(1) = theta_sum(j);
w(1) = w_sum(j);
a(1) = a_sum(j);
b(1) = b_sum(j);

% Update Potential Energy
h = (1 – cos(theta(1)));
% h can be b(1)
%h = length*(1 – b(1));
h = length * h;
Potential(j) = mass*g*h;

% Update Kinetic Energy
v = length * w(1);
Kinetic(j) = 0.5*mass*v*v;

% Update Total Energy
Total_Energy(j) = Potential(j) + Kinetic(j);


% Plot Results

% Plot Theta Results
plot(1:Frames, theta_sum,’b-‘);
title(‘Parker-Sochacki Implemetation of Simple Pendulum’);
xlabel(‘Frame Number’);

% Plot w Results
% plot(1:Frames, w_sum,’b-‘);
% title(‘Parker-Sochacki Implemetation of Simple Pendulum’);
% xlabel(‘Frame Number’);
% ylabel(‘w’);

% Plot Potential Energy Results
% plot(1:Frames, Potential,’b-‘);
% title(‘Potential Energy of Simple Pendulum’);
% xlabel(‘Frame Number’);
% ylabel(‘Potential Energy’);

% Plot Kinetic Energy Results
% plot(1:Frames, Kinetic,’b-‘);
% title(‘Kinetic Energy of Simple Pendulum’);
% xlabel(‘Frame Number’);
% ylabel(‘Kinetic Energy’);

% Plot All Potential and Kinetic Energies
% plot(1:Frames, Potential,’r’, 1:Frames, Kinetic, ‘b’);
% legend(‘Potential Energy’, ‘Kinetic Energy’);
% title(‘Potential Energy and Kinetic Energy of a Simple Pendulum’);
% xlabel(‘Frame Number’);
% ylabel(‘Energy’);

% Plot All Energies
% plot(1:Frames, Potential,’r’, 1:Frames, Kinetic, ‘b’, 1:Frames, Total_Energy, ‘o-‘);
% legend(‘Potential Energy’, ‘Kinetic Energy’, ‘Total Energy’);
% title(‘Potential Energy and Kinetic Energy of a Simple Pendulum’);
% xlabel(‘Frame Number’);
% ylabel(‘Energy’);

% Plot Total Energy
% plot(1:Frames, Total_Energy, ‘o-‘);
% title(‘Total Energy of Simple Pendulum’);
% xlabel(‘Frame Number’);
% ylabel(‘Total Energy’);


By just looking at the plot, the pendulum swings back and forth without appearing to gain or lose energy-not significantly anyway. The time step chosen was 1/24 mimicking the frame rate for a movie. The order of the Maclaurin series is 19. This means that at an order of 19, the pendulum is swinging for over 20 seconds in a film without significant changes(we would have to look at several decimal points to see where the values are changing).

Now I will lower the order of the Maclaurin series for theta, w, and the auxiliary variables. When I lowered the order to a 2nd order, the simulation quickly blew up. Even when I lowered the time step from 1/24 to 1/240. you can see the sine wave getting larger meaning the total energy is gaining meaning that the swing of the pendulum is getting larger:


The timing is also off because this should be over the same amount of time as the first simulation. This means that in the first simulation the pendulum came back to its original position 9 times and in this version the simulation came back to its original position(and overshot it) 10 times.

To better see what is going on, we need to look at plots in terms of energy.Because Total Energy should be conserved, Potential Energy should hit its maximum value when Kinetic Energy is at its minimum value. Also, Kinetic Energy should hit its maximum value when Potential Energy is at its minimum value. In a perfect setting, neither should gain over time. Here is a plot of Potential Energy and Kinetic Energy over the first setting, where the time step is 1/24 and the order is 19:PE and KE_good.PNG

As opposed to second order and a time step of 1/240:

PE and KE_Bad.PNG

The Kinetic Energy is growing more rapidly. If this simulation were to continue, it would soon blow up.

Keep in mind that this is a second order polynomial at a time step 1/10 of the 19th order method and the energy is still gaining!

Let’s look at a plot of Total Energy and see the difference. Here is the first set of conditions where the order is 19 and the time step is 1/24:


In theory, we would love to see a perfect line and the total energy is getting bigger as time goes on. But look at the range of this simulation. The total energy stays between values of 10.4999999999998 and 10.5000000000002. This is likely coming from round off error but this shows that Parker-Sochacki is conserving the total energy of the simple pendulum when it is at a higher order.

Look at the second setting, where the order is 2 and the time step is lowered to 1/240:


The Total Energy will soon triple.

These were extremely different values.

Starting at 3rd order, Here is a list of the minimum and maximum energy values going out to only 4 digits(I went to many more at 19th order).

Order    Total Energy Minimum    Total Energy Maximum   Difference

3                     10.4998                                    11.1517                             0.6519

4                     10.4266                                    10.50                                0.0734

5                      10.5004                                    10.4977                           0.0027

6                      10.5004                                    10.50                                0.0004

7                       10.5000                                    10.5000                          0.0000

Unless I go to more digits, we won’t see a pattern beyond a 7th order Maclaurin series.

On an interesting note, Orders 3, 6, and 7 were gaining energy. Orders 4 and 5 were losing energy. But as we increased the order, the total energy was better conserved.

This pattern was similar if I changed initial conditions. When I increased the length of the pendulum or changed the value of the mass of the bob or changed the value of gravity, total energy was better conserved as order increased. The same held true under different initial conditions for theta and w.

As another test, I increased the time step to 10/24 and raised the order of the Maclaurin series to 99:


The minimum Total Energy was 10.4998 and the maximum Total Energy was 10.5008. While eventually the time step will be too large where Total Energy cannot be conserved, for this example, a much larger time step could be taken and have the Total Energy be conserved.

This shows a lot of of promise. It is possible to have the Parker-Sochacki method conserve energy at larger time steps for the Simple Pendulum when the order is raised. While this was only applied to the Simple Pendulum, it shows promise that energy will be conserved in other physical simulations that conserve energy.