코드는 교수님 강의자료에 첨부된 코드를 참고했습니다.

close ALL

clear

clc

%% Mass

m=[4 4 4];

k=[4 4 4];

W=[1 0];

M=diag(m);

x0=[1;0;0];

v0=[0;0;0];

t=0:0.01:10;

%% Global Stiffness Matrix

elem=size(k,2);

nodes=elem+1;

K=zeros(nodes);

% Connectivity Matrix

for i=1:elem

 con(i,1:2)=[i i+1];

end

for i=1:elem

 ke=k(i)*[ 1 -1;

 -1 1];

 node1=con(i,1);

 node2=con(i,2);

 K(node1,node1)=K(node1,node1)+ke(1,1);

 K(node1,node2)=K(node1,node2)+ke(1,2);

 K(node2,node1)=K(node2,node1)+ke(2,1);

 K(node2,node2)=K(node2,node2)+ke(2,2);

end

if W(1)==1 && W(2)==0

 K=K(2:end,2:end);

elseif W(1)==0 && W(2)==1

 K=K(1:end-1,1:end-1);

elseif W(1)==1 && W(2)==1

 K=K(2:end-1,2:end-1);

else

 K=K;

end

%% Calcaultions

L=sqrt(M);

Kmn=inv(L)*K*inv(L);

[v,lamda]=eig(Kmn);

P=v; % check P`*P=I

S=inv(L)*P;

r0=inv(S)*x0;

rdot0=inv(S)*v0;

w=diag(sqrt(lamda));

for i=1:size(m,2)

 r(i,:)=sqrt(r0(i)^2+rdot0(i)^2/w(i)^2)*cos(w(i)*t); % atan값이 pi/2이므로 cos(w(i)*t)로 입력

end

x=S*r;

%% Plot

figure(1)

plot(t,x(1,:),`linewidth`,2)

hold on

plot(t,x(2,:),`linewidth`,1)

plot(t,x(3,:),`linewidth`,1)

legend_entries={`x1(t)`,`x2(t)`,`x3(t)`};

legend(legend_entries,``,`northeast`);

grid on

xlabel(`Time (s)`)

ylabel(`Displacement (m)`)

title(`Three-Degree-of-Freedom Model (Undamped)`)  

49번째 줄 r(i,:)=sqrt(r0(i)^2+rdot0(i)^2/w(i)^2)*sin(w(i)*t+atan(w(i)*r0(i)/rdot0(i))); 에서

r(i,:)=sqrt(r0(i)^2+rdot0(i)^2/w(i)^2)*cos(w(i)*t); 로 바꿔서 진행했습니다.

그 결과 산출된 그래프도 같이 첨부했습니다.

데스크탑으로 돌려도 원하는 그래프가 나오지 않아서 sqrt(r0(i)^2+rdot0(i)^2/w(i)^2) 를 r0(i)로 바꿨더니 올바른 답이 나왔습니다.

그래서 이 부분에 오류가 있다고 생각했습니다.

제가 잘못 생각한 부분이 있으면 어디서 생긴 문제인지 알려주세요.

감사합니다.