matlab求微分方程组数值解出现问题

Question

我已经反复检查该方程无误 但是在ode45求解这个微分方程组的时候解出来的值全是NAN，（以下是流程图）论文里说龙格库塔能解啊。

clear all
%输入参数
gam = 0.0728;
rou = 1000;
g = 9.8;
b = 0.006;
a = (gam/(rou*g))^(1/2);
theta = 90*pi/180;
V_0 = 25 *10^(-6);
p = 2*a/b;
c = 10^(-11);
D = 2*((3*V_0/(4*pi))^(1/3));
Bo =(rou*g*D^2)/gam;
%在theta小于135°以及Bo大于1的时候求解
if theta<135*pi/180&&Bo>1
    %这个是要求解的方程
dxdt = @(t,y) y(1) * cos(t) / (y(1)*y(2) + p*y(1) - sin(t));
dzdt = @(t,y) y(2) * sin(t) / (y(1)*y(2) + p*y(1) - sin(t));

tspan = [0 theta];
x0 = 0;
z0 = 0;
%设置收敛条件 用V_final的一个范围，或者b已经降到了0下
flag = true;
  while flag == true
[t, y] = ode45(@(t, y) [dxdt(t, y); dzdt(t, y)], [0, theta], [x0; z0]);

final_state = y(end,:);
final_time = t(end);
V_final = abs(pi*final_state(1)*(final_state(1) * final_state(2) + p * final_state(1) - 2 * sin(final_time)));
    if (V_final > (0.9*V_0)/(a^3) && V_final < (1.1*V_0)/(a^3))||b<0
   flag = false;
    else
 %如果这样的b无法满足上面的收敛条件，那么就微调b
  b = b -0.00001;
  flag = true;
    end

  end
%输出
X = a*y(:,1);
Y = a*final_state(2) -a*y(:,2);
plot (X,Y)
else 
%下面应该这些没问题
    R = (3*V_0/pi*((cos(theta)^3)-3*cos(theta)+2))^(1/3);
    center = [0,R];
    alpha = linspace(0, theta, 100);
    X= a* R * sin(alpha)  ;
    Y = a*center(2) -a*R * cos(alpha); 
    axis equal;
plot (X,Y)

end

但是一旦给把初始值设成很小的量，就能解出结果，但是结果显然不正确（谁家水滴是分子量级的高度啊）愁死。。
以防万一我把方程贴上

Answer 1

2301_79402751 晚上好🌙🌙🌙
本答案参考ChatGPT-3.5

根据你的描述，你在使用ode45求解微分方程组时出现了问题。以下是一些可能的解决方案：

1. 检查初始值的合理性：确保初始值适当，不要设置为非常小的值，因为这可能导致数值不稳定。可以尝试使用与问题相匹配的初始值进行解决。

2. 检查方程的正确性：仔细检查方程的推导过程和变量的定义，并确保没有错误。可能会出现符号错误、变量定义错误或方程形式错误等问题。

3. 检查函数的定义：确保函数dxdt和dzdt被正确定义，输入参数的顺序和数量与方程中的一致。

4. 调整ode45的参数：可以尝试调整ode45的参数，例如减小相对误差容限（RelTol）或绝对误差容限（AbsTol），以提高数值稳定性。

5. 检查方程中的除零错误：在方程中存在除以零的情况，可能会导致NaN的结果。请确保方程中分母不为零，或在计算之前添加相应的判断。

6. 对比论文中的求解方法：仔细阅读论文中关于求解微分方程组的方法，并对比已有的代码实现，确保求解方法正确。

7. 尝试其他求解方法：如果ode45仍然无法解决问题，可以尝试其他求解方法，如ode15s、ode23、ode23s等。

希望以上解决方案对你有帮助。如果问题仍然存在，请提供更多细节或代码，以便进一步帮助你解决问题。