警告: 在 t=4.465718e-05 处失败。在时间 t 处，若不将步长降至允许的最小值(1.084202e-19)以下，积分公差要求无法满足

Question

matlab中如何解决:警告: 在 t=4.465718e-05 处失败。在时间 t 处，若不将步长降至允许的最小值(1.084202e-19)以下，积分公差要求无法满足的问题
   这是我的代码
m = 0.015; 
Fa = 70; 
g = 9.81; 
c = 0.2; 
Gr = 0.004; 
n = 1500; 
alpha = 40 * pi / 180; 
Z = 7; Ri = 2.945; Re = 3.06; Dw = 5.678; Di = 12; De = 32; 
Bd = Ri + Re - Dw + Gr; 
z0=0.001;
zi=z0-Bd*sin(alpha);
ze=z0+Bd*sin(alpha);

% 材料参数
nu1 = 0.3; nu2 = 0.3; 
epsilon1 = 2.03 * 10^(5); epsilon2 = 2.03 * 10^(5); 
rhoi = -(2 / Di + 1 / Ri); 
rhoe = -(2 / De + 1 / Re); 

E = 1.0003 + 0.5968 * (rhoe / rhoi); 
F = 1.5277 + 0.6203 * log(rhoi / rhoe); 
k = 1.0339 * (rhoi / rhoe)^0.636; 
epsilon = 2 * epsilon1 * epsilon2 / (epsilon2 * (1 - nu1^2) + epsilon1 * (1 - nu2^2)); 
Ki = pi^2 * k^2 * epsilon^2 * E / (4.5 * F^3 * rhoi);
Ke = pi^2 * k^2 * epsilon^2 * E / (4.5 * F^3 * rhoe);
Ki = abs(Ki);
Ke = abs(Ke);
K = 1 / (Ki^(-2/3) + Ke^(-2/3));

% 时间设置
tspan = [0, 0.002]; 
y0 = [0.001; 0; 0.001; 0; 0.001; 0]; 

omegae = n * 2 * pi / 60; 
omegac = omegae/2;
omegai= 0;

[t, y] = ode45(@(t, y) zhendong(t, y, m, c, K, Bd, g,omegac, Fa, omegai, Z, Re, Ri, zi, ze), tspan, y0);

function dydt = zhendong(t, y, m, c, K, g, Bd, Fa, omegac, omegai, Z, Re, Ri, zi,  ze)
    dydt = zeros(6, 1); 
    Fx = 0;
    Fy = 0;
    Fz = 0;


     Niw = 7;
     New = 7;
     BOi0 = 0.01; BOip = 0.2;
     BOe0 = 0.01; BOep = 0.3;
     lambdai = pi * Ri / Niw;
     lambdae = pi * Re / New;
   
    for j = 1:Z
        alphaem = 2 * pi * (j - 1) / Z + omegac * t;
         alphaim = 2 * pi * (j - 1) / Z + (omegac - omegai) * t; 
    
         Li = Ri * alphaim;
         Le = Re * alphaem;
        
         % 内圈和外圈波纹度
         BOi = BOi0 + BOip * sin(2 * pi * Li / lambdai);
         BOe = BOe0 + BOep * sin(2 * pi * Le / lambdae);

        hx = (Re - Ri) * cos(alphaem) + y(1);
        hy = (Re - Ri) * cos(alphaem) + y(3);
        hz = ze-zi+y(5);

        Bd1 = sqrt(hx^2 + hy^2 + hz^2);
        
        delta = Bd1 - Bd + BOi + BOe; 
        phi_j = hx / Bd1;
        psi_j = hy / Bd1;
        chi_j = hz / Bd1;
                
        cosphi_j=cos(phi_j);
        cospsi_j=cos(psi_j);
        coschi_j=cos(chi_j);
        
        Fx = Fx + K * abs(delta)^1.5 * cosphi_j;
        Fy = Fy + K * abs(delta)^1.5 * cospsi_j;
        Fz = Fz + K * abs(delta)^1.5 * coschi_j;
    end

    dydt(1) = y(2); 
    dydt(2) = -(c * y(2) + Fx) / m; 
    dydt(3) = y(4); % 
    dydt(4) = (m * g - c * y(4) - Fy) / m; 
    dydt(5) = y(6); 
    dydt(6) = (Fa - c * y(6) - Fz) / m; 
end

Answer 1

本答案参考 多次询问 ChatGPT-3.5 后选择的最佳回答，若是满意，还请采纳，谢谢

1. 问题分析
   • 这个警告表明在使用ode45求解常微分方程（ODE）时，在特定的时间t = 4.465718e - 05处，为了满足积分公差要求，需要将步长降低到非常小的值（1.084202e - 19）以下，但这是不被允许的。这可能是由于问题的刚性（stiffness）或者函数zhendong中的某些计算导致数值不稳定。
2. 可能的解决方案
   • 调整公差参数
     • 在ode45中，可以尝试调整相对公差（RelTol）和绝对公差（AbsTol）。默认情况下，RelTol为1e - 3，AbsTol为1e - 6。例如，可以将相对公差设置得更宽松一些：

       options = odeset('RelTol',1e - 2);
       [t, y]=ode45(@(t, y) zhendong(t, y, m, c, K, Bd, g,omegac, Fa, omegai, Z, Re, Ri, zi,  ze), tspan, y0,options);
       

   • 检查函数zhendong中的计算
     • 数值稳定性问题
       • 在zhendong函数中，计算Fx、Fy和Fz的循环部分可能存在数值不稳定的情况。例如，abs(delta)^1.5这种计算可能在delta接近0时导致数值问题。可以考虑添加一些条件判断来避免这种情况，比如当delta接近0时，采用线性近似或者其他更稳定的计算方式。
     • 参数取值和单位一致性
       • 检查所有的参数（如m、c、K等）的取值是否合理，并且确保在整个计算过程中单位是一致的。例如，如果某些参数的单位不匹配，可能会导致计算结果在数值上出现异常，进而影响到ode45的求解。
   • 尝试其他求解器
     • 如果问题是刚性的，ode45可能不是最佳选择。可以尝试使用ode15s（专门用于求解刚性问题的求解器）：

       options = odeset('RelTol',1e - 3,'AbsTol',1e - 6);
       [t, y]=ode15s(@(t, y) zhendong(t, y, m, c, K, Bd, g,omegac, Fa, omegai, Z, Re, Ri, zi,  ze), tspan, y0,options);