The program

We transform the second order differential equation to a system of two first order differential equations

To resolve this system, we use the 4th order Runge-Kutta method

The complete matlab program is shown above :

fichier instab.m

x(1)=input('enter the initial position x(1)= ');

y(1)=input('enter the initial velocity y(1)= ');

f =input('enter the forcing amplitude f= ');

w =input('enter the forcing frequency w= ');

b =input('enter the dampening constant b = ');

N =input('enter the nombre of iterations N = ');

z(1)=f;

t(1)=0;

h=pi/(50*w); % Time step which allows us to do 100 times

% steps per drive cycle(one complete period of

%the forcing term)

k=1;

% Variables to extract values to plot the poincar section

xx(1)=x(1);

yy(1)=y(1);

%Mean loop to resolve with the Runge-Kutta method

for i = 1 : N

x1=h*fx(y(i));

y1=h*fy(x(i),y(i),t(i),f,b,w);

x2=h*fx( y(i)+.5*y1 );
y2=h*fy(x(i)+.5*x1, y(i)+.5*y1, t(i)+.5*h,f,b,w);

x3=h*fx( y(i)+.5*y2);

y3=h*fy(x(i)+.5*x2, y(i)+.5*y2, t(i)+.5*h,f,b,w);

x4=h*fx( y(i)+y3);

y4=h*fy(x(i)+x3, y(i)+y3, t(i)+h,f,b,w);

x(i+1) = x(i) + (1/6)*(x1 + 2*x2 + 2*x3 +x4);

y(i+1) = y(i) + (1/6)*(y1 + 2*y2 + 2*y3 + y4);

t(i+1) = t(i)+h;

z(i+1) = f*cos(w*t(i+1));

%Condition to extract only the values of x and y corresponding

%to a value of z(t)=f*cos(w*t)=f so w*t is a multiple of 2*pi

% so t(i)*w=m*2*pi=(i-1)*h -> i-1=k*2*pi/(h*w)-> i=100*m+1 so

%the remainder of the division of i per 100 must be equal to 1

if rem(i+1,100)==1

k=k+1;

xx(k)=x(i+1);

yy(k)=y(i+1);

end
end

%position vs time plot

figure(1);

plot(t,x);

xlabel('time');

ylabel('position');

title('position vs time');

%velocity vs time plot

figure(2);

plot(t,y);

xlabel('time');

ylabel('velocity');

title('velocity vs time');

%The phase space plot(2d)

figure(3);

plot(x,y);

xlabel('position');

ylabel('velocity');

title('phase space');

%The true phase space plot(3d)with the forcing term

figure(4);

plot3(x,y,z);

xlabel('position');

ylabel('velocity');

zlabel('forcing');

title('three dimensionnal phase space')

%The poincare section plot

figure(5);

plot(xx,yy,'go');%plot only a green circle at the data point

%without tracing a line between them

title('poincaré section');

xlabel('position');

ylabel('velocity');

fichier fx.m

function fonct1= fx(y)
fonct1=y;

Fichier fy.m

function fonct2= fy(x,y,t,f,b,w)
fonct2= x - x^3 - b*y + f*cos(w*t);

This program permit us to make different plots:

-Position vs time x(t)

-Velocity vs time y(t)

-The phase space velocity vs position y=f(x)

-The true phase space which is three dimensionnal by adding the forcing term z(t)=f*cos(wt)

-The poincaré section