我正在网上阅读这篇文章,其中提到使用“if语句”和“abs()”函数会在MATLAB的变步长ODE求解器(如ODE45)中产生负面影响.根据OP,它可以显着影响时间步长(需要太低的时间步长),并且当微分方程最终被积分时给出差的结果.我想知道这是否属实,如果是,为什么.此外,如果不诉诸固定步骤求解器,如何减轻这个问题.我在下面给出了一个示例代码,我的意思是：

function [Z,Y] = sauters(We,Re,rhos,nu_G,Uinj,Dinj,theta,ts,SMDs0,Uzs0,...

Uts0,Vzs0,zspan,K)

Y0 = [SMDs0;Uzs0;Uts0;Vzs0]; %Initial Conditions

options = odeset('RelTol',1e-7,'AbsTol',1e-7); %Tolerance Levels

[Z,Y] = ode45(@func,zspan,Y0,options);

function DY = func(z,y)

DY = zeros(4,1);

%Calculate Local Droplet Reynolds Numbers

Rez = y(1)*abs(y(2)-y(4))*Dinj*Uinj/nu_G;

Ret = y(1)*abs(y(3))*Dinj*Uinj/nu_G;

%Calculate Droplet Drag Coefficient

Cdz = dragcof(Rez);

Cdt = dragcof(Ret);

%Calculate Total Relative Velocity and Droplet Reynolds Number

Utot = sqrt((y(2)-y(4))^2 + y(3)^2);

Red = y(1)*abs(Utot)*Dinj*Uinj/nu_G;

%Calculate Derivatives

%SMD

if(Red > 1)

DY(1) = -(We/8)*rhos*y(1)*(Utot*Utot/y(2))*(Cdz*(y(2)-y(4)) + ...

Cdt*y(3)) + (We/6)*y(1)*y(1)*(y(2)*DY(2) + y(3)*DY(3)) + ...

(We/Re)*K*(Red^0.5)*Utot*Utot/y(2);

elseif(Red < 1)

DY(1) = -(We/8)*rhos*y(1)*(Utot*Utot/y(2))*(Cdz*(y(2)-y(4)) + ...

Cdt*y(3)) + (We/6)*y(1)*y(1)*(y(2)*DY(2) + y(3)*DY(3)) + ...

(We/Re)*K*(Red)*Utot*Utot/y(2);

end

%Axial Droplet Velocity

DY(2) = -(3/4)*rhos*(Cdz/y(1))*Utot*(1 - y(4)/y(2));

%Tangential Droplet Velocity

DY(3) = -(3/4)*rhos*(Cdt/y(1))*Utot*(y(3)/y(2));

%Axial Gas Velocity

DY(4) = (3/8)*((ts - ts^2)/(z^2))*(cos(theta)/(tan(theta)^2))*...

(Cdz/y(1))*(Utot/y(4))*(1 - y(4)/y(2)) - y(4)/z;

end

end

功能“dragcof”由以下内容给出：

function Cd = dragcof(Re)

if(Re <= 0.01)

Cd = (0.1875) + (24.0/Re);

elseif(Re > 0.01 && Re <= 260.0)

Cd = (24.0/Re)*(1.0 + 0.1315*Re^(0.32 - 0.05*log10(Re)));

else

Cd = (24.0/Re)*(1.0 + 0.1935*Re^0.6305);

end

end