%clear all
m=1;
ro=1.2;
J0=1.;
k=4;
c=2;
a=-0.01;
S=1.;
L=1;
V=3.;
T=50;
NPT=2000;
dt=T/NPT;
z0=0;
dz0=.1;
alpha0=0.03;
dalpha0=0.;
beta0=0.;
%
temp=linspace(1,T+1,NPT)-1;
beta=linspace(1,2,NPT)*0.;
z=linspace(1,2,NPT)*0.;
alpha=linspace(1,2,NPT)*0.;
%========
% Construction de M K et F
%========
M=[m, a*cos(alpha0);a*cos(alpha0),J0];
K=[k ,-ro*S*V^2*dcz0(alpha0)/2;0 c-ro*S*L*dcm0(alpha0)/2];
R=linsolve(M,K);
%
F0=[ro*S*V^2*cz0(alpha0),ro*S*L*V^2*cm0(alpha0)/2];
F=linsolve(M,F0');
%==================
% donnees initiales
%==================
z(1)=z0;
z(2)=z0(1)+dt*dz0(1);
alpha(1)=alpha0;
alpha(2)=alpha(1)+dt*dalpha0;
X1=[z(1),alpha(1)]';
X2=[z(2),alpha(2)]';
%=============
% resolution du systeme edo
%===========
for i=3:NPT
X3=2*X2-X1-dt^2*R*X2+dt^2*F;
z(i)=X3(1);
alpha(i)=X3(2);
X1=X2;
X2=X3;
end
%=======
figure
plot(temp,z,temp,alpha)