我在网上阅读了这篇文章,其中有人提到使用“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
最佳答案
这是因为使用 if
语句、模运算 ( abs()
) 或单位阶跃函数、狄拉克增量等计算的导数会在解或其导数的值中引入不连续性,从而导致扭结、跳跃、拐点等。
这意味着 ODE 的解决方案在相关时间完全改变了行为。可变步长积分器将做的是
因此,问题点附近会出现许多失败的步骤和步长的减小,从而对整体积分时间产生负面影响。
然而,可变步长积分器将继续产生良好的结果。恒定步长积分器不是解决此类问题的好方法,因为它们首先无法检测到此类问题(没有错误估计)。
您可以做的只是将问题分成多个部分。如果您事先知道将在什么时间点发生变化,您只需为每个间隔开始新的积分,每次使用前一次积分的输出作为下一次积分的初始值。
如果您事先不知道问题出在哪里,您可以在 Matlab 的 ODE 求解器中使用这个非常好的功能,称为事件函数(参见 the documentation )。您让 Matlab 的求解器之一检测事件(导数中的符号变化、
if
语句中的条件变化或其他),并在检测到此类事件时终止积分。然后开始新的积分,从上次开始,并使用前一次积分的初始条件,和以前一样,直到达到最终时间。这样总体执行时间仍然会有轻微的损失,因为 Matlab 会尝试准确检测事件的位置。但是,就执行时间和结果的准确性而言,它仍然比盲目运行集成要好得多。
关于Matlab:可变步长 ODE 求解器中的 if 语句和 abs() 函数,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/16582841/