please see this code and tel me where is the problem beacause i don't find it

hold on

%Nouveau programme le 21/02/2014

%%%%%%%%%%%%%Données physique du probleme%%%%%%%%%%%%%%%%%%%%%%%%%%

M=30;                                   % masse de la poutre KG/ML??
EI=0.20*10^7;                            % Rigidité flexionnelle de la poutre N/m^2
ksi=0.4;
k=0.04*10^8;                             % rigidité du sol N/m^2
v=27.77;                                %vitesse m/s de chargement
n1=5;
c=0.05;

%nombre des modes

l=2*(3*pi)/2*(4*EI/k)^0.25;               %m;longueur elastique de la poutre en m

%%%%%%%%%%%%%%%% Les conditions initiales %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

v0=0.001;                               %m/s;vitesse initiale
T0=0.0001;                              %m déplacement initiale
t1=l/v;                                 %temps de parcours de la longueuer elastique  
ww1=0;
Q=1;                              %chargement en Newton

n=23;                                 % Nombre de points d'analyse

f=-21:1:21;                          %creation du vecteu englobant toute la longueur élastique 

beta=1;

for m=1:5

    wn(m)=((((m^4)*(pi^4)*EI)/(M*(l^4)))+(k/M))^0.5;   % Pulsation propre du systeme
    wd(m)=(wn(m)^2-beta^2)^0.5;                      % pseudo pulsation
    omega(m)=m*pi*v/l; 
for i=1:n
    f(i)=i;
    xx(i)=l/(n-1);
        xx1(i)=xx(i)*(i-1);

phi(i,m)=sin(m*pi*(xx1(i))/l);

for j=1:n
  tt(j)=l/v/(n-1);
           tt1(j)=tt(j)*(j-1);
          x1(j)=v*tt1(j) 
A1=2*Q/(M*l);

A2(m)=1/(wd(m)^4+2*((beta^2-wd(m)^2)*omega(m)^2)+omega(m)^4); %(((w_d(m))^4)+2*((((ksi*w_d)^2)+((w_n(m)^2)))*(mu(m)^2))+mu(m)^4);

A3(m)=2*beta*omega(m);

A4(m)=(omega(m)/wd(m))*(((beta)^2)-(wd(m)^2)+omega(m)^2);                              

A5(m)=-A3(m);

A6(m)=wd(m)^2-omega(m)^2;

T11(j,m)=exp(-beta*(tt1(j)));

T1(j,m)=A1*A2(m)*(T11(j,m)*(A3(m)*cos(wd(m)*tt1(j))+A4(m)*sin(wd(m)*tt1(j)))+A5(m)*cos(omega(m)*tt1(j))+A6(m)*sin(omega(m)*tt1(j)));

ww1(i,j,m)=phi(i,m)*T1(j,m);    

end
end
end

for i=1:n;
  for j=1:n;
      ww1(i,j)=0;
      for m=1:n1;
ww1(i,j)=ww1(i,j)+ww1(i,j,m);

%ww1(i,j) = sum(ww1(i,j,:));

        end

    end
end

for i=1:n
  www(i)=ww1(i,j)
  j=5
end
%plot(www,'g')

for j=1:n
  www1(j)=ww1(i,j)
  i=12
end
plot(www1,'r')    
%plot(ww1(10,:))

[Walter Roberson]

What difference is there between what you expect and what shows up when you run the program?

[sihem chaib]

thank you to be Interest; Yes the curve is breaken and i have to plot it fonction of x the beam lenght not fonction the number of element 'i' and

Answer 1

[Youssef Khmou]

The program is functioning well , here is the result :

[sihem chaib]

yes but i dn't have the result i want the cuve is breaken

[Youssef Khmou]

you mean the result interpretation? ok the code is technically correct, the problem then is in the constants, verify the constants /or change : c, ksi and the velocity v .