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 

Я не могу получить код для запуска, я получаю следующее сообщение об ошибке: решение Не найдено ни одного.

Я пробовал устанавливать допуски всего 1e-20, но это ничего не меняло.

Answer 1

Как ваша система настроена, на самом деле удобно ее строить и наблюдать за ее поведением. Я векторизовал вашу функцию и построил f(x) = MLNEfun(x)-x, где выход MLNE(x) - newA. Вы заинтересованы в фиксированной точке вашей системы.

Я наблюдаю это:

Существует особенность, и корень скрещивание, при А ~ 3800. Мы можем использовать fzero, так как это квадратные скобки корень решатель, и дать ему очень жесткие границы на решении fzero(@(x)MLNEfun(x)-x, [3824,3825]), который производит 3.8243e+03. Это пара порядков величин от вашего исходного предположения. Нет решения для вашей системы около ~ 3e5.

Update В спешке, мне не удалось увеличить на участке, который показывает другой (хорошо себя) корень в 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