MATLAB中ode45报错：需返回列向量问题求助

错误原因分析 #

1. 返回向量类型不符合要求：ode45求解器要求状态方程的匿名函数必须返回列向量，但原代码中用逗号分隔元素，生成的是行向量，直接触发odearguments报错。
2. 强迫力维度不匹配：预先生成的Forcing是与time等长的数组，而ode45会在任意时间点进行求解计算，直接调用Forcing/m会导致维度不匹配的问题。

具体修复步骤 #

1. 调整匿名函数返回格式 #

将匿名函数中各状态导数的分隔符从逗号,改为分号;，确保返回的是列向量，符合ode45的输入要求。

2. 实时计算强迫力 #

删除预先生成的Forcing数组，改为在匿名函数内根据当前时间t直接计算强迫力：F*cos(omega*t)/m，彻底避免维度不匹配问题。

3. 优化代码可读性（可选） #

将状态方程拆分到多行（用...换行），避免一行代码过长，方便后续调试和维护。

修改后的完整代码 #

%% 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数组，改为实时计算

%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); ...
        F*cos(omega*t)/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]);
    % 可添加结果保存或绘图代码，避免循环覆盖t,y
    % 示例：save(['result_case_',num2str(j),'.mat'],'t','y');
end

额外说明 #

原代码中每次循环会覆盖变量t和y，如果需要保留所有工况的计算结果，建议在循环内将结果存入数组或单独保存为文件。

内容的提问来源于stack exchange，提问作者Quentin Anderson-Watson