Programmation d'un vecteur

%% Données
E = 50000000000; %Béton  (en GPa)
rho = 1500; %Béton 
e = 4; %épaisseur de la pont (m)
l = 2; %largeur de la poutre supportant la chaussée (m)
S = e*l; %(en m2)
m_v = 1200; %masse moyenne des véhicules(kg)
g = 9.81; %accélération de la pesanteur
I = l*e^3/12; %moment quadratique
v = 50; %vitesse moyenne des véhicules (km/h)

%% Paramètres spatio-temporels
tf = 3600;  %Temps final de l'étude (en s) 
Npast = 10;    %Nombre de pas de temps
dt = tf / Npast; %Pas de temps
T = linspace(0,tf,Npast);

L = 270;   %Longueur de la portion étudiée (en m)
Npasx = 10;  %Nombre de pas d'espace
dx = L / Npasx;   %Pas d'espace
X = linspace(0,L,Npasx);

%% Constantes diverses
a = rho*S*dx^4;
b = E*I*dt^2;
c = 6*b-2*a;
d = (dx^4)*(dt^2);

%%Conditions initiales
mu = 0;
gamma = 1;

%% Création de la matrice A
A = [];
for i=2:Npasx-1
    A(i,i+1) = 4*b;
    A(i+1,i) = 4*b;
    A(i+2,i) = -b;
    A(i,i+2) = -b;
    A(i,i) = -c;
end
A(1,1) = -c-b;
A(Npasx,Npast)= -c+4*b;
A(Npasx,Npast-1)= 3*b;

%% Création du vecteur h
h = [Npasx];
for j=1:dt:tf
    if (dx*i-v*dt*j)< dt && (dx*i-v*dt*j)>-dt
        h(j)=0
    else h(j)=1
    end
end  
%% Création du vecteur p
p = [];
for j=1:dt:Npast
    for i=1:dx:Npasx
    p(j) = rho*g*S + m_v*g*h*(dx*i-v*dt*j)
    end
end

%% Création de la matrice w
w = zeros(Npasx,Npast);
w(:,0) = mu;
w(:,1) = gamma*dx + mu;

for k=1:Npast
    w(2:Npasx-1,k) =(1/a)*( A*w(2:Npasx-1,k-1) - a*w(2:Npasx-1,k-2) - d*p'); 
end

mesh(T,X,w)
grid on
grid minor
o = title ('Méthode numérique')
p = xlabel('Temps (en s)')
q = ylabel('Position dans la poutre (en m)')
r = zlabel('Déplacements verticaux(en m)')

o.FontSize = 18
p.FontSize = 15
q.FontSize = 15
r.FontSize = 15