我正在尝试使用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 ，但这并没有改变任何东西 .
1 回答
您的系统设置方式，绘制它并观察其行为实际上很方便 . 我向你的函数进行了矢量化并绘制了
f(x) = MLNEfun(x)-x，其中MLNE(x)的输出是newA. 实际上，您对系统的固定点感兴趣 .我观察到的是：
在A~3800处有一个奇点和一个根交叉 . 我们可以使用
fzero，因为它是一个括号中的根求解器，并且在解决方案fzero(@(x)MLNEfun(x)-x, [3824,3825])上给它非常紧的界限，它产生了3.8243e+03. 从您的开始猜测开始，这是几个数量级 . 您的系统在~3e5附近没有解决方案 .Update 在我的匆忙中，我没有放大图，它显示了另一个（表现良好的）根在1.3294e 04.由你决定哪个是物理上有意义的根 . 我在下面说的一切仍然适用 . 在您感兴趣的解决方案附近开始猜测 .
In response to comment 由于你想为
K的不同值执行此操作，那么你最好的办法就是坚持fzero，只要你要求一个变量而不是fsolve. 这背后的原因是fsolve使用牛顿方法的变体，这些方法没有括号，并且很难在像这样的奇点找到解决方案 . 另一方面，fzero使用布伦特的方法，该方法保证在括号内找到根（如果存在） . 它在奇点附近表现得更好 .在执行布伦特方法之前，MATLAB的
fzero的实现也会搜索包围间隔 . 因此，如果您提供足够接近根的单个起始猜测，它应该为您找到它 . 以下fzero的样本输出：