# Calculating Lyapunov exponent (LE) for pendulum using ellipsoid growth - code yields negative LEs

I was redirected here from physics stack exchange where hopefully my question is more appropriate. Per my advisor, I have read the textbook Chaos, an introduction to dynamical systems by Alligood, Sauer, and Yorke. (side note, I have really enjoyed this book). In chapter 5, the numerical calculation of Lyapunov exponents (LE) is given where you track the growth of the ellipsoid using the Jacobian of the system and Gram-Schmidt orthogonalization. Example 5.6 asks one to determine the LEs for the Henon map with parameters $$a = 1.4$$ and $$b = 0.3$$ using this method, noting one should expect to receive LEs of $$h_1 = 0.42$$ and $$h_2 = -1.62$$, respectively, for the dimensions $$x$$ and $$y$$. I have successfully done this in Matlab so I know the code works.

My issue is that when I try to apply this to the forced, damped pendulum, I output negative LEs when input parameters should yield chaotic behavior. I have used the equations from My Physics Lab Chaotic Pendulum so I know what parameters should yield this chaotic behavior (assuming this site evaluated it accurately). It is solved using ODE45. The initial conditions I input are random.

The inputs for the Jacobian I use are those that result from the time-$$2\pi$$ map. When that failed, I tried using data for every iteration instead. Bust.

According to the site and my code, I should be seeing a positive LE, but both I get are negative. Any advice? Thank you!

Here is the code. The henon map is also included, but commented out. Currently, it is set to using the time-$$2\pi$$ map data, but one could easily switch it out for every iteration instead:

% henon
% chaos an intro to dyn, yorke, pg 201

% housekeeping
clear
clc
format compact
close all

% %% Henon data
% % parameters
% a = 1.4;
% b = 0.3;
% % iterations
% N = 10^3;
% %allocating space for variables
% x = zeros(N,2);
% L = zeros(N,2);
% % orthogonal basis for R2 map
% W = eye(2);
%
% % IC
% x(1,:) = rand(2,1);
%
% %iterate the map
% for i =1:N
% x(i+1,1) = a-x(i,1)^2+b*x(i,2);
% x(i+1,2) = x(i,1);
% end
%
% %plot those points
% plot(x(:,1),x(:,2),'.');

%% Pendulum Data

% site: https://www.myphysicslab.com/pendulum/chaotic-pendulum-en.html

% driving 'frequency', but really w
k = 2/3;
% Poincare plot will sample at this driving frequency
% T = 1/f = 1/w/(2*pi) --> T = 2*pi/w
T_f = 2*pi/k;
%natural freq = sqrt(g/R)
g = 1;
m = 1;
b = 0.5;
W = eye(2);
% F_d values for certain pendulum behaviors:
% single: 0.9, 1.35
% double: 1.07, 1.45
% chaotic: 1.15, 1.50

F_d = 1.5;

t_begin = 0;
t_end = 10000;

%theta0 = [0 0]; %IC
theta0 = rand(1,2);
[ts, thetas] = ode45(@(t,theta) pend_damp_force(t,theta,k,g,Radius,m,b,F_d), [t_begin t_end], theta0);

% to get time-2pi map
time_points = [0:2*pi/k:ts(end)];
time_points = dsearchn(ts, time_points');

pos_poin = wrapToPi(thetas(:,1));
vel_poin = thetas(:,2);
pos_2pi_map = pos_poin(time_points);
vel_2pi_map = vel_poin(time_points);
N = length(pos_2pi_map);
%% LE calc

%matrix times circle = ellipse, for every point
for j = 1:N
% Henon
%     J = [-2*x(j,1),b; % jacobian
%         1,0];
% Pendulum
J = [0 1;
% https://www.math.ucla.edu/~yanovsky/Teaching/Math151B/handouts/GramSchmidt.pdf
[W,R] = qr(J*W); % QR is same thing as gran-schmidt...A = QR where Q is orthogonal matrix (Q^TQ=I) and R is upper triangular matrix
L(j,:) = diag(R);
end

%calc LE and display
LE = mean(log(abs(L)),1) % the one is to do it along each column, if it were 2 it would do it along each row


The function pend_damp_force:

function dthetadt = pend_damp_force(t,theta,k,g,R,m,b,F_d)

dthetadt(2) = -g/R*sin(theta(1))+ (-b*theta(2)+F_d*cos(k*t))/(m*R^2);  % g(pos,vel)

end

• Have you tested out that small differences in the initial conditions lead to eventually very different solutions? Commented Jan 27, 2022 at 11:00
• Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
– Community Bot
Commented Jan 27, 2022 at 16:57
• @LutzLehmann, yes I have confirmed that the system when input with parameters that the site states cause chaotic behavior is sensitive to initial conditions and its time-2pi map results in bounded behavior. Also, although the bot stated this, I am unsure how to make my question more clarified? Unless it is about the ellipse/gram-schmidt method, but I assumed anyone answering would have to already have knowledge of that technique. Commented Jan 27, 2022 at 20:20
• I haven't checked, but maybe a sign error for the Jacobian? Commented Sep 15, 2022 at 10:19

Good morning, good afternoon or good night :-). I made two programs in octave on the subject of forced pendulum, as described in my physics lab (website, Chaotic Pendulum: https://www.myphysicslab.com/pendulum/chaotic-pendulum-en.html). In this first code I use the dynamic equations of the forced pendulum to plot the phase space (angle x angular velocity). In case you have any doubts, I put Taylor's book on classical mechanics below in bibliographies, where I also did the research.

    %==============================================================================
%{
FEDERAL UNIVERSITY OF GOIÁS (UFG)
INSTITUTE OF PHYSICS (IF-UFG)

SUBJECT                   : COMPUTATIONAL PHYSICS
STUDENT                   : CARLOS EDUARDO DA SILVA LIMA
DATE                      : 11/01/2021
MODIFIED FOR THE LAST TIME: 11/07/2021
THEME                     : FORCED PENDULUM
LANGUAGE/SOFTWARE         : OCTAVE
DATA ENTRY                : ENTER THE DATA IN THE SECTIONS MARKED WITH THE SYMBO
L, IN THIS SECTION, THE USER CAN CHANGE THE VALUES.
PURPOSE                   : IN THIS CODE, WE CALCULATE VIA 4 ORDER RUNGE-KUTTA
THE PARAMETERS T, X1(T) AND X2(T) AND WE PLOT THE GR
APHICS X1(T) VS X2(T), THE SYSTEM PHASE SPACE, LIMIT
ADO IN A RANGE FROM -PI TO PI. FIGURE 1, REPRESENTS
IN THE TEXT ARTICLE.
%}
% ==============================================================================
% ==============================================================================
clc
clear
disp('FORCED PENDULUM');
% ==============================================================================
% PARAMETERS, FORCED PENDULUM AND 4TH ORDER RUNGE-KUTTA.
% - IN THIS PART, THE USER CAN CHANGE THE PARAMETERS.
global q  = 0.05;     % Q = K/M.
global b  = 1.02;     % F0/(M*L).
global w0 = 2/3;      % ANGLE FREQUENCY RAD/S OF EXTERNAL DRIVING FORCE.
t_0(1) = 0;          % INITIAL TIME INSTATE.
x1_0(1) = (pi/180)*60;% INITIAL ANGLE.
x2_0(1) = 0;          % INITIAL ANGLE SPEED.
h = 0.01;             % STEP, OR SIZE OF SUB INTERVAL.
N = 500000;           % NUMBER OF ITERATIONS .
% ==============================================================================
% ORDINARY DIFFERENTIAL EQUATIONS IN THE FIRST ORDER - 2X2 EDO SYSTEM.

% ODE 1.
function f1 = f1(t,x1,x2)
f1 = x2;
endfunction

% ODE 2.
function f2 = f2(t,x1,x2)
global q;
global b;
global w0;
f2 = -q*x2 - sin(x1) + b*cos(w0*t);
endfunction
% ==============================================================================
% ITERATION VIA 4 ORDER RUNGE-KUTTA.
t(1) = t_0; % INITIAL TIME INSTATE.
x1(1) = x1_0; % INITIAL ANGLE.
x2(1) = x2_0; % INITIAL ANGLE SPEED.

for i = 1:1:(N-1)
K11 = h*f1(t(i), x1(i), x2(i));
K12 = h*f2(t(i), x1(i), x2(i));
K21 = h*f1(t(i) + (h/2), x1(i) + (K11/2), x2(i) + (K12/2));
K22 = h*f2(t(i) + (h/2), x1(i) + (K11/2), x2(i) + (K12/2));
K31 = h*f1(t(i) + (h/2), x1(i) + (K21/2), x2(i) + (K22/2));
K32 = h*f2(t(i) + (h/2), x1(i) + (K21/2), x2(i) + (K22/2));
K41 = h*f1(t(i) + h, x1(i) + K31, x2(i) + K32);
K42 = h*f2(t(i) + h, x1(i) + K31, x2(i) + K32);

x1(i+1) = x1(i) + ((K11 + 2*(K21 + K31) + K41)/6);
x2(i+1) = x2(i) + ((K12 + 2*(K22 + K32) + K42)/6);
t(i+1) =  t(i) + h;
endfor

% WE LIMIT THE VALUES OF X1(T) AND X2(T) BETWEEN THE RANGE [-PI,PI].
for i = 1:1:N
if x1(i)>-pi && x1(i)<pi
teta(i) = x1(i);
w(i)    = x2(i);
endif
endfor

% DATA OUTPUT ON SCREEN.
for i = 1:1:N
% printf("t(%d) = %f || teta(%d) = %f || w(%d) = %f\n",i,t(i),i,x1(i),i,x2(i));
endfor
% ==============================================================================
% GRAPHICS.

% FIGURE 1
figure(1)
plot(teta, w,'b.');
title('Espaço de fase para pêndulo forçado.')
h=get(gcf, "currentaxes");
set(h, "fontsize", 20, "linewidth", 1.5);
grid()

%===============================================================================
%                   BIBLIOGRAPHIES AND RESEARCH SITES
% SITE: https://www.myphysicslab.com/pendulum/chaotic-pendulum-en.html
% BOOK: Classical Mechanics: John R. Taylor


In this second part, I reuse the first code inserted above to calculate the Lyapunov exponents (Graph 1) and the sensitivity of the system regarding the imposed conditions (Graph 2). In this code, I simulate two pendulums with the same attributes (same mass, same length of the rods etc), but with slight differences in the initial conditions. In this simulation I chose to keep the initial velocities of both forced pendulums equal (Both pendulums start from rest with certain initial angles), but in the second pendulum I make a small increase in the radian eps value (eps is any small value, where eps = 0.02). I will leave the text in Brazilian South American Portuguese, I hope I have helped you as well as the community. Note: Please, if you disagree with something, let me know.

% SEGUNDO PÊNDULO
y1(1) = x1(1) + 0.02;
y2(1) = x2(1);


Code 2

% ==============================================================================
%{
INSTITUTO DE FÍSICA (IF-UFG)

DICIPLINA                 : FÍSICA COMPUTACIONAL
ALUNO                     : CARLOS EDUARSDO DA SILVA LIMA
DATA                      : 01/11/2021
LINGUAGEM/SOFTWARE        : OCTAVE
(*), NESTA SEÇÃO, O USUÉRIO PODE ALTERAR OS VALORES.
OBJETIVO                  : NESTE CÓDIGO, CÁLCULAMOS VIA RUNGE-KUTTA DE 4 ORDEM
OS PARÃMETROS T, X1(T) E X2(T) E PLOTAMOS OS GRÁFICO
S T VS LOG, ALÉM DE T VS X1 E T VS X2.
%}
% ==============================================================================
% ==============================================================================
clc
clear
% ==============================================================================
% PARÃMETROS DO PÊNDULO FORÇADO E RUNGE-KUTTA DE 4ORDEM
% (*) - NESTA PARTE, O USUÁRIO PODE ALTERAR OS PARÂMETROS.
global q  = 0.05; % Q = K/M.
global b  = 1.02; % F0/ML.
global w0 = 2/3;  % FREQUÊNCIA ÂNGULAR RAD/S DA FORÇA MOTRIZ EXTERNA.
t_0(1) = 0;           % TEMPO INICIAL
x1_0(1) = (pi/180)*60; % ÂNGULO INICIAL
x2_0(1) = 0;           % VELÔCIDADE INICIAL
h  = 0.01;  % PASSO (TAMANHO DAS SUB-DIVISÕES DO INTERVALO)
N  = 10000; % QUANTIDADE DE ITERAÇÕES
% ==============================================================================
% EQUAÇÕES DIFERENCIAIS ORDINÁRIAS NA PRIMEIRA ORDEM - SISTEMA 2X2 EDO.
% EDO 1
function f1 = f1(t,x1,x2)
f1 = x2;
endfunction
% EDO 2
function f2 = f2(t,x1,x2)
global q;
global b;
global w0;
f2 = -q*x2 - sin(x1) + b*cos(w0*t);
endfunction
% ==============================================================================
% ITERAÇÃO VIA RUNGE-KUTTA DE 4 ORDEM
% PRIMEIRO PÊNDULO
t(1) =  t_0;% TEMPO INICIAL
x1(1) = x1_0;% ÂNGULO INICIAL
% SEGUNDO PÊNDULO
y1(1) = x1(1) + 0.02;
y2(1) = x2(1);
for i = 1:1:(N-1)

% PRIMEIRO PÊNDULO
K11 = h*f1(t(i), x1(i), x2(i));
K12 = h*f2(t(i), x1(i), x2(i));
K21 = h*f1(t(i) + (h/2), x1(i) + (K11/2), x2(i) + (K12/2));
K22 = h*f2(t(i) + (h/2), x1(i) + (K11/2), x2(i) + (K12/2));
K31 = h*f1(t(i) + (h/2), x1(i) + (K21/2), x2(i) + (K22/2));
K32 = h*f2(t(i) + (h/2), x1(i) + (K21/2), x2(i) + (K22/2));
K41 = h*f1(t(i) + h, x1(i) + K31, x2(i) + K32);
K42 = h*f2(t(i) + h, x1(i) + K31, x2(i) + K32);

% SEGUNDO PÊNDULO
Kl11 = h*f1(t(i),y1(i),y2(i));
Kl12 = h*f2(t(i),y1(i),y2(i));
Kl21 = h*f1(t(i)+(h/2), y1(i)+(Kl11/2), y2(i)+(Kl12/2));
Kl22 = h*f2(t(i)+(h/2), y1(i)+(Kl11/2), y2(i)+(Kl12/2));
Kl31 = h*f1(t(i)+(h/2), y1(i)+(Kl21/2), y2(i)+(Kl22/2));
Kl32 = h*f2(t(i)+(h/2), y1(i)+(Kl21/2), y2(i)+(Kl22/2));
Kl41 = h*f1(t(i)+h, y1(i) + Kl31, y2(i) + Kl32);
Kl42 = h*f2(t(i)+h, y1(i) + Kl31, y2(i) + Kl32);

% PRIMEIRO PÊNDULO
x1(i+1) = x1(i) + ((K11 + 2*(K21 + K31) + K41)/6);
x2(i+1) = x2(i) + ((K12 + 2*(K22 + K32) + K42)/6);

% SEGUNDO PÊNDULO
y1(i+1) = y1(i) + ((Kl11 + 2*(Kl21 + Kl31) + Kl41)/6);
y2(i+1) = y2(i) + ((Kl12 + 2*(Kl22 + Kl32) + Kl42)/6);

% TEMPO, IGUAL PARA OS DOIS PÊNDULOS FORÇADOS.
t(i+1) = t(i) + h;
endfor

% LOG - EXPOENTE DE LYAPUNOV, LOG|X1-Y1|
for i = 1:1:N
Log(i) = log10(abs(x1(i)-y1(i)));
endfor
% ==============================================================================
% GRAFICOS

% FIGURA 1
figure(1)
plot(t, Log,'b.');
title('Expoente de Lyapunov.')
xlabel('Tempo(s)'), ylabel('Log|Ângulo|')
h=get(gcf, "currentaxes");
set(h, "fontsize", 20, "linewidth", 1.5);
grid()

% FIGURA 2
figure(2)
plot(t,x1,'b.',t,y1,'r.');
title('Pêndulo 1 (Azul), Pêndulo 2 (Vermelho).')