fsolve/fzero: 未找到解决方案，显示正常

Question

我正在尝试使用 fsolve 或 fzero 执行以下算法：

K5=8.37e-2
P=1

Choose an A
    S2=(4*K5/A)^(2/3)
    S6=3*S2
    S8=4*S2

    SO2 = (5*P)/149 - (101*S2)/149 - (293*S6)/149 - (389*S8)/149
    H2O = (40*P)/149 + (556*S2)/447 + (636*S6)/149 + (2584*S8)/447
    H2S = 2*SO2

    newA = (H2O)^2/(SO2)^3

Repeat until newA=oldA

主要解决的是K5=1/4 * A * S2^3/2。 S2就是从这里算出来的

所以这是我在 Matlab 中所做的：

function MultipleNLEexample
clear, clc, format short g, format compact

Aguess = 300000; % initial guess

options = optimoptions('fsolve','Display','iter','TolFun',[1e-9],'TolX',[1e-9]); % Option to display output
xsolv=fsolve(@MNLEfun,Aguess,options); 

[~,ans]=MNLEfun(xsolv)

%- - - - - - - - - - - - - - - - - - - - - -
function varargout = MNLEfun(A); 
    K5 = 8.37e-2;

    S2 = (4*K5/A)^(2/3);
    S6 = 3*S2;
    S8 = 4*S2;

    P=1; %atm

    SO2 = (5*P)/149 - (101*S2)/149 - (293*S6)/149 - (389*S8)/149;
    H2O = (40*P)/149 + (556*S2)/447 + (636*S6)/149 + (2584*S8)/447;
    newA=H2O^2/SO2^3;
    fx=1/4*newA*S2^(3/2)-K5;

    varargout{1} = fx;
    if nargout>1
        H2S = 2*SO2;
        varargout{2} = ((2*S2+6*S6+8*S8)/(2*S2+6*S6+8*S8+H2S+SO2)*100);
    end

我无法将我的代码获取到运行，我收到以下错误：未找到解决方案。

fsolve 已停止，因为根据梯度测量问题看起来很规则，但函数值的向量不接近于零，由函数公差的选定值。

我尝试将公差设置为低至 1e-20，但这并没有改变任何东西。

Answer 1

你的系统设置方式，实际上很方便绘制它并观察它的行为。我矢量化了你的函数并绘制了 f(x) = MLNEfun(x)-x，其中 MLNE(x) 的输出是 newA。实际上，您对系统的 不动点 感兴趣。

我观察到的是：

在 A ~ 3800 处有一个奇异点和一个根交叉。我们可以使用 fzero 因为它是一个带括号的根求解器，并在解 fzero(@(x)MLNEfun(x)-x, [3824,3825]) 上给它非常严格的边界产生 3.8243e+03。这与您的初始猜测相差几个数量级。 ~3e5附近没有你系统的解决方案。

更新匆忙中，我未能放大绘图，它显示了 1.3294e+04 处的另一个（表现良好的）根。由您决定哪个是物理上有意义的。我在下面所说的一切仍然适用。只需在您感兴趣的解决方案附近开始猜测即可。

回复评论 由于您想对 K 的不同值执行此操作，因此最好的选择是坚持使用 fzero，只要您求解的是一个变量，而不是 fsolve。这背后的原因是 fsolve 使用牛顿法的变体，这些变体没有括号，并且会像这样在奇异点上很难找到解决方案。 fzero 另一方面，使用 Brent 的方法，保证在括号内找到根（如果存在）。它在奇异点附近也表现得更好。

MATLAB 的 fzero 实现也在执行布伦特方法之前搜索包围间隔。因此，如果您提供一个足够接近根的起始猜测，它应该会为您找到它。下面 fzero 的示例输出：

fzero(@(x)MLNEfun(x)-x, 3000, optimset('display', 'iter'))

Search for an interval around 3000 containing a sign change:
 Func-count    a          f(a)             b          f(b)        Procedure
    1            3000       -616789          3000       -616789   initial interval
    3         2915.15       -433170       3084.85       -905801   search
    5            2880       -377057          3120  -1.07362e+06   search
    7         2830.29       -311972       3169.71  -1.38274e+06   search
    9            2760       -241524          3240  -2.03722e+06   search
   11         2660.59       -171701       3339.41  -3.80346e+06   search
   13            2520       -109658          3480  -1.16164e+07   search
   15         2321.18      -61340.4       3678.82   -1.7387e+08   search
   17            2040      -29142.6          3960   2.52373e+08   search

Search for a zero in the interval [2040, 3960]:
 Func-count    x          f(x)             Procedure
   17            2040      -29142.6        initial
   18         2040.22      -29158.9        interpolation
   19         3000.11       -617085        bisection
   20         3480.06  -1.16224e+07        bisection
   21            3960   2.52373e+08        bisection
   22         3720.03  -4.83826e+08        interpolation
....
   87         3824.32  -5.46204e+48        bisection
   88         3824.32   1.03576e+50        bisection
   89         3824.32   1.03576e+50        interpolation

Current point x may be near a singular point. The interval [2040, 3960] 
reduced to the requested tolerance and the function changes sign in the interval,
but f(x) increased in magnitude as the interval reduced.

ans =

   3.8243e+03

fsolve/fzero: 未找到解决方案，显示正常

fsolve/fzero: No solution found, appears regular

matlab

mathematical-optimization

numerical-methods

solver

nonlinear-functions