The program

To solve numerically the differential systems in the pendulum equation, we just need to discretize it. I had first chosen an easy Euler (explicit) scheme for each equations. It appears that this scheme is unstable for a damping factor equal to zero. With the recommendation of A.Merle, I changed the scheme, and the second equation is solved with an implicit scheme. The results seem to be very good.

The coding resulting is simple and it is like coding a serie.

This program is written in matlab language and is suitable for each case developed in this project. In the beginning, the variable type allow us to choose the type of pendulum we want to simulate. Then, several parameters can be modified, such as the damping factor and the excitation.

clear;

type=input('Type de pendule a modeliser (1=normal, 2=avec forcage): ');

nbcas=input('Nombre de cas a tester: ');

n=input('Nb de pas de tps simules: ');

Dt=0.01;
w=9.81;

for j=1:nbcas

u(1,j)=input('Angle de depart: ');
v(1,j)=input('Vitesse de depart: ');

t(1)=0;

if type == 1

l=input('Valeur de l ammortissement: ');

for i=1:n
u(i+1,j)=v(i,j)*Dt + u(i,j);
v(i+1,j)= 1/(1+Dt*l) * (-Dt*w*sin(u(i+1,j)) + v(i,j));
t(i+1)=t(i)+Dt;
end

end

if type == 2

m=input('Valeur du forcage: ');
k=input('Valeur k de l ammortissement (de forme k*teta^2): ');

for i=1:n
u(i+1,j)=v(i,j)*Dt + u(i,j);
v(i+1,j)= 1/(1+Dt*(k*u(i+1,j)*u(i+1,j)-m)) * (-Dt*w*sin(u(i+1,j)) + v(i,j));

t(i+1)=t(i)+Dt;
end

end

end

plot(u,v);
figure;
plot(t,u);

The outputs of the program are the phase portrait plot(u,v) and the evolution of the angle u during the time plot(t,u).