matlabode45

MATLAB error in ode45, says that it must return a column vector

Error using odearguments @(T,Z)[Z(2),K_RUBBER/M_PLYWOODZ(3)+C_RUBBER/M_PLYWOODZ(4)-(C_RUBBER+C(J))/M_PLYWOODZ(2)-(K_RUBBER+K(J))/M_PLYWOODZ(1),Z(4),FORCING/M+C_RUBBER/MZ(2)+K_RUBBER/MZ(1)-C_RUBBER/MZ(4)-K_RUBBER/MZ(3)] must return a column vector.

%% Project Parameters
%{
    Isolation mats use 5% of the floor area (A_mat=.05A)
%}

clc,clear

time=linspace(0,3,5001);            %Time in seconds
A = (8*12*.0254)^2;                 %Area of floor
A_mat = .05*A;                      %Area of mat
m_plywood = (1200)*(A)*(3*.0254);   %Mass of plywood
m = 550/9.81;                       %Mass of runner
m0 = 1e-3*min([m_plywood m]);       %Mass separating isolation mats

%Parameters for Harmonic Forcing
T =72/60;       %Period of forcing from runner               
f = 1/T;        %Frequency of forcing from runner
omega = f*2*pi; %Frequency in radians of forcing from runner
F = -m*9.81;    %Force runner applies to floor
Forcing = F*cos(omega*time);

%Parameters for rubber mat
m_rubber = (1800)*A*(1*.0254);                      %Rubber mat mass
k_rubber = (12e6)*A/(1*.0254);                      %Rubber mat stiffness
zeta_rubber = .08;                                  %Rubber mat damping ratio
c_rubber = zeta_rubber*2*sqrt(k_rubber*m_rubber);   %Rubber mat damping coefficient

%Parameters for SPX
m_spx = (1200)*A_mat*(1*.0254);           %SPX mat mass
k_spx = (1e6)*A/(1*.0254);                %SPX mat stiffness
zeta_spx = .08;                           %SPX mat damping ratio
c_spx = zeta_spx*2*sqrt(k_spx*m_spx);     %SPX mat damping coefficient

%Parameters for DMP
m_dmp = (1500)*A_mat*(1*.0254);           %DMP mat mass
k_dmp = (10e6)*A/(1*.0254);               %DMP mat stiffness
zeta_dmp = .22;                           %DMP mat damping ratio
c_dmp = zeta_dmp*2*sqrt(k_dmp*m_dmp);     %DMP mat damping coefficient

%Parameters for FRX
m_frx = (1800)*A_mat*(1*.0254);           %FRX mat mass
k_frx = (8e6)*A/(1*.0254);                %FRX mat stiffness
zeta_frx = .15;                           %FRX mat damping ratio
c_frx = zeta_frx*2*sqrt(k_frx*m_frx);     %FRX mat damping coefficient

%Parameters for BPX
m_bpx = (1200)*A_mat*(1*.0254);           %BPX mat mass
k_bpx = (6e6)*A/(1*.0254);                %BPX mat stiffness
zeta_bpx = .06;                           %BPX mat damping ratio
c_bpx = zeta_bpx*2*sqrt(k_bpx*m_bpx);     %BPX mat damping coefficient

k = [k_bpx, k_frx, k_dmp, k_spx]';        %Setting stiffness vector to loop through
c = [c_bpx, c_frx, c_dmp, c_spx]';        %Setting damping vector to loop through

%% Dynamic System Modeling, Scenario 1: 1 inch of ioslation mats
clc
for j = 1:length(k)
    [t,y]=ode45(@(t,z)[z(2),k_rubber/m_plywood*z(3)+c_rubber/m_plywood*z(4)-(c_rubber+c(j))/m_plywood*z(2)-(k_rubber+k(j))/m_plywood*z(1),z(4),Forcing/m+c_rubber/m*z(2)+k_rubber/m*z(1)-c_rubber/m*z(4)-k_rubber/m*z(3)],time,[0;0;0;0]);
end
Solution

  • This function handle returns a row vector since you use commas as separators:

    @(t,z)[z(2),k_rubber/m_plywood*z(3)+c_rubber/m_plywood*z(4)-(c_rubber+c(j))/m_plywood*z(2)-(k_rubber+k(j))/m_plywood*z(1),z(4),Forcing/m+c_rubber/m*z(2)+k_rubber/m*z(1)-c_rubber/m*z(4)-k_rubber/m*z(3)]

    Use semi-colons instead to make this return a column vector:

    @(t,z)[z(2);k_rubber/m_plywood*z(3)+c_rubber/m_plywood*z(4)-(c_rubber+c(j))/m_plywood*z(2)-(k_rubber+k(j))/m_plywood*z(1);z(4);Forcing/m+c_rubber/m*z(2)+k_rubber/m*z(1)-c_rubber/m*z(4)-k_rubber/m*z(3)]

    For code readability, it would probably be better to define this function handle outside of your ode45( ) call.