MATLAB: 2 Degree of Freedom system with ODE 45

2 dofode 45

I am trying to solve a 2 DOF system using ODE 45, and plot the displacement and velocity response. I believe I am very close but my velocity graph isn't showing up as expected. The problem may be in my initial condition matrix or my EOM function file. The initial conditions are supposed to be x1=.2, x2=.1, v1=v2=0. I've messed around with the placement of the IC's in the matrix to try and get the right response.
function [Xdot] =EOM(tspan,X,k1,k2,k3,c1,c2,c3,m1,m2,F0,w) 
Xdot =zeros(4,1);
x1=X(1);        
x1dot=X(2);      
Xdot(1)=x1dot; 
x2=X(3);
x2dot=X(4);
Xdot(2)=x2dot;
Xdot(2,1)= (-((k1+k2)*x1)/m1)+((k2*x2)/m1)-(((c1+c2)*x1dot)/m1)+((c2*x2dot)/m1)+((F0*cos(w*tspan))/m1); 
Xdot(4,1)= (-((k2+k3)/m2)*x2)+((k2/m2)*x1)-(((c2+c3)*x2dot)/m2)+((c2*x1dot)/m1);
end
clc
clear all
close all
m1=4;
m2=9;
k1=10;
k2=10;
k3=10;
c1=2;
c2=2;
c3=2;
F0=.5;
w=20;
EOM0=@(tspan,X)EOM(tspan,X,k1,k2,k3,c1,c2,c3,m1,m2,F0,w); 
X0=[.2 0 .1 0];  
tspan=[0 10];
[T,X] = ode45(EOM0,tspan,X0);  
figure(1)
plot(T,X(:,2))
hold on
plot(T,X(:,4))
title('Displacement with Damping and Harmonic Force')
xlabel('Time, \it t, \rm (s)')
ylabel('Displacement, \it x, \rm (m)')
legend('X1','X2')
figure(2)
[T,X] = ode45(EOM0,tspan,X0);  
plot(T,X(:,1))
hold on
plot(T,X(:,3))

Best Answer

Try
function [Xdot] =EOM(tspan,X,k1,k2,k3,c1,c2,c3,m1,m2,F0,w) 
Xdot =zeros(4,1);
x1=X(1);        
x1dot=X(2);     
x2=X(3);
x2dot=X(4); 
Xdot(1,1)=x1dot; 
Xdot(2,1)= (-((k1+k2)*x1)/m1)+((k2*x2)/m1)-(((c1+c2)*x1dot)/m1)+((c2*x2dot)/m1)+((F0*cos(w*tspan))/m1); 
Xdot(3,1)=x2dot;
Xdot(4,1)= (-((k2+k3)/m2)*x2)+((k2/m2)*x1)-(((c2+c3)*x2dot)/m2)+((c2*x1dot)/m1);
end
Best wishes
Torsten.
Related Question