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，因为它是一个括号根求解器，并给它很生产3.8243e+03的解决方案fzero(@(x)MLNEfun(x)-x, [3824,3825])的紧密边界。这是你开始猜测的几个量级。 ~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