MATLAB嵌套for循环调用ode45时无法运行求助

问题描述 #

我在求解多组合弹簧-质量-阻尼系统的动力学问题，单层for循环下代码运行完全正常，但嵌套for循环无法运行，且无任何报错信息，ode45函数无法完成求解，恳请提供解决建议。

原MATLAB代码 #

%% 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 seperating isolation mats

m=158.757; % kg of person and weight
% impact velocity is 4.89 m/s for 200 lb object dropped from 48 inches
% by impulse momentum theorem, we have ( 90.7 kg )*( 4.89 m/s) = ( M=158.757 kg ) * ( inital velocity )
% solving for initial velocity, we get 2.7937 m/s
xdot0=90.7*4.89/158.757; % 2.7937 m/s initial velocity of floor

%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,y1]=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); ...
        c_rubber/m*z(2)+k_rubber/m*z(1)-c_rubber/m*z(4)-k_rubber/m*z(3)],time,[0;xdot0;0;0]);
end


%% Dynamic System Modeling, Scenario 2: 2 inches of isolation mats

clc
for j = 1:length(k)
   for jj = 1:length(k)
        [t,y2]=ode45(@(t,z)[z(2); ...
            c(jj)/m0*z(4)+k(jj)/m0*z(3)-(c(jj)+c(j))/m0*z(2)-(k(jj)+k(j))/m0*z(2); ...
            z(4); ...
            c_rubber/m_plywood*z(6)+k_rubber/m_plywood*z(5)-(c_rubber+c(jj))/m_plywood*z(4)-(k_rubber+k(jj))/m_plywood*z(3); ...
            z(6); ...
            c_rubber/m*z(4)+k_rubber/m*z(3)-c_rubber/m*z(6)-k_rubber/m*z(5); ...
            ],[0:15],[0;xdot0;0;0;0;0]);
    end
end

%% Dynamic System Modeling, Scenario 3: 3 inches of isolation mats
clc
for j = 1:length(k)
    for jj = 1:length(k)
        for jjj = 1:length(k)
            [t,y3]=ode45(@(t,z)[z(2); ...
                c(jj)/m0*z(4)+k(jj)/m0*z(3)-(c(jj)+c(j))/m0*z(2)-(k(jj)+k(j))/m0*z(1); ...
                z(4); ...
                c(jjj)/m0*z(6)+k(jjj)/m0*z(5)-(c(jjj)+c(jj))/m0*z(4)-(k(jjj)+k(jj))/m0*z(3); ...
                z(6); ...
                c_rubber/m_plywood*z(8)+k_rubber/m_plywood*z(7)-(c_rubber+c(jjj))/m_plywood*z(6)-(k_rubber+k(jjj))/m_plywood; ...
                z(8); ...
                c_rubber/m*z(6)+k_rubber/m*z(5)-c_rubber/m*z(8)-k_rubber/m*z(7); ...
                ],time,[0;xdot0;0;0;0;0;0;0]);
        end
    end
end

问题排查与修复建议 #

• 修正状态方程笔误

  • Scenario 2中第二个方程的刚度项错误：-(k(jj)+k(j))/m0*z(2) 需改为 -(k(jj)+k(j))/m0*z(1)，否则系统动力学方程完全偏离物理模型，引发数值求解发散。
  • Scenario 3中第五个方程缺失位移项：-(k_rubber+k(jjj))/m_plywood 需改为 -(k_rubber+k(jjj))/m_plywood*z(5)，否则方程维度不匹配，触发隐性数值异常。

• 统一时间输入参数
  Scenario 2中使用[0:15]作为ode45的时间输入，与Scenario 1、3的time参数区间、输出点密度不一致，建议统一替换为time，保证求解条件一致。

• 更换刚性求解器
  代码中m0为极小值（1e-3倍最小质量），导致系统方程出现极大系数，使系统成为刚性系统。ode45是针对非刚性系统的求解器，对刚性系统求解效率极低甚至无法完成，建议替换为ode15s或ode23t等刚性求解器。

• 修复结果存储逻辑
  嵌套循环中每次迭代都会覆盖t,y2、t,y3变量，最终仅保留最后一组结果，且反复赋值可能引发内存波动。建议使用cell数组或三维数组存储所有组合结果，示例：

  % Scenario 2初始化存储
  y2 = cell(length(k),length(k));
  t2 = cell(length(k),length(k));
  for j = 1:length(k)
     for jj = 1:length(k)
          [t2{j,jj},y2{j,jj}]=ode15s(...); % 替换为刚性求解器
      end
  end
  

• 移除冗余清屏命令
  嵌套循环内的clc会清空控制台输出，无法观察求解过程中的隐性警告或中间信息，建议仅保留开头的clc，或直接移除以便排查问题。

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