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

 

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)로 바꿨더니 올바른 답이 나왔습니다.

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

 

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

 

감사합니다.​