MATLAB: How to solve and plot the equation of motion with ODE solver? (Keep getting errors)

odevibration

I am trying to solve and plot my equation of motion for a vibration system. The original equation of motion is: ((m*l^2/6)+M*l^2)*theta_dd + c*theta_d + (E*b*h^3/6*l^3)*theta = F*cos(omega*t).
Here is the code I have written so far:
function ligo_simulation
% Define initial conditions x(0) & x'(0)
x0 = [0 1];
% Time interval
t = [0 60];
% Solve equation of motion
[T, X] = ode45(@EOM, t, x0);
plot(T, X(:,1), T, X(:,2));
% EOM function & variables
function dx = EOM(t, x)
E = 187.6E9;  % Modulus of Elasticity for stainless steel (Pa)
b = 0.216;    % Base (m)
h = 13.87E-3; % Height (m)
l = 0.429;    % Length (m)
M = 1224;     % Hanging mass (kg)
zeta = 0.05;  % Damping ratio
f = 5;        % Desired frequency (Hz)
% Smaller mass (m*l^2/6) is negligible compared to hanging mass because it
% is magnitudes less
m = M*l^2;               % Mass (kg)
k = ((E*b*h^3)/(6*l^3)); % Spring constant
c = zeta*2*(sqrt(k*m));  % Damping constant
omega_n = sqrt(k/m);     % Natural frequency
omega = 2*pi*f;          % Frequency
F = [0 100];
term1 = F*cos(omega*t);
term2 = k*x(1);
term3 = c*x(2);
term4 = m;
dx = zeros(2,1);
dx(1) = x(2);
dx = [x(2) ((term1 - term2 - term3)/(term4))]; % Equation of motion
fprintf('The natural frequency is %0.2f rad/s \n', omega_n)
end
end
And here is the error message I keep getting:
Error using odearguments (line 93)
LIGO_SIMULATION/EOM must return a column vector.
Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in ligo_simulation (line 11)
[T, X] = ode45(@EOM, t, x0);
I've tried changing my variable around and I can't figure out why it says I am not returning a column vector. If I change things from [0 60] to [0; 60], I get a concatenating error. Eventually I want to vary my variables with a for loop and graph a bunch of values, but I need to make this work first. Any help is much appreciated. Thanks!

Best Answer

Hi Eric,
There are two things that keep this from running. First, you initialize dx as a column vector, but the statement
dx = [x(2) ((term1 - term2 - term3)/(term4))]; % Equation of motion
makes a row vector. It's missing the semicolon
 dx = [x(2); ......]
as you found. That would fix things up except that the second part involving term1, term2 etc. is not a scalar. That's because F = [0 100], a 1x2 row vector, and that property carries through the rest of the algebra, so you end up trying to concatenate a 1x2 row vector below the scalar x(2). I'm not sure what your intent was for F, but changing to F = 100 gives code that runs and looks pretty good.
Related Question