前言

非线性方程组，顾名思义就是未知数的幂除了不是1，其他都有可能！线性方程组其实只是非常小的一类，非线性方程组才是大类！也正因此非线性方程组包含各种各样的方程形式，所以它的解和对应的求解方法不可能像线性方程组那样完美，即都是局部收敛的。先给出一个直观的非线性方程组例子：

$\begin{cases} x_{1}^{2} - 10*x_{1} + x_{2}^{2} + 8 = 0 & \\ x_{1}*x_{2}^{2} + x_{1} - 10*x_{2} + 8 = 0 & \end{cases}$

个人对两个问题的理解：

1. 非线性方程组如果有解，一般都有很多解！如何理解：

把方程组的解看成是各个函数图像的交点的。我们知道非线性方程组的各个函数就都是复杂曲线/面，甚至是高纬空间里的复杂东西；线性方程组的各个函数就是最简单的直线/面！各个复杂函数图像间的相交机会很多，并且只要相交，就是多个交点(因为交线、交面里有无数的交点)，也就是有多个解~

可以想象，非线性方程组有多解是很平常的一件事~ 对于复杂的非线性函数没解才不正常！
可以想象，这些解是等价的！没有说是等级更高，谁等级低一些。都是解！因为：只要是解，它就只满足一个条件：让方程组中的各个方程=0。所以无法用什么评判标准(比如范数)来说哪个解的等级高一些或者效果更好一些。
注意：这里的解等价和欠定线性方程组通解中的唯一极小范数解不一样！可以想象二者的区别：非线性方程组中的解都是实打实存在的；而欠定线性方程组中除了特解，其他通解中的解说存在也行，说不存在那就是因为方程条件(个数)都不够！这些是啥都行的通解和非线性方程组中实打实存在的解肯定不能比！

这样的话各个非线性方程组的局部收敛性就可以理解，即：空间中有很多解时，我每次只能找一个，那我找谁？找离我出发点最近的那个解呗。所以不同的出发点，就有可能找到不同的解，这就是局部收敛性。

2. 为什么不能用简单的替换先变成线性方程组求解，然后再解每一个非线性方程？

意思是：每个方程中把所有和 $x_1$ 有关的用一个变量代替，所有 $x_2$ 有关的用一个变量代替，即方程1中用： $w_1 = x_{1}^{2} - 10*x_{1}$ ， $w_2 = x_{2}^{2} + 8$
但是很明显方程2的第一项 $x_{1}*x_{2}^{2}$ 两个变量相乘，没法用变量代替~
并且，即使在方程2中能代替，那么就会有 $w_3$ 和 $w_4$ ，这样总未知数变成4个而方程只有2个，还是解不了。
所以，非线性方程组不可能用简单的线性变量代换来解。

本文介绍最常用、最好用的非线性方程组解法，包括牛顿法和拟牛顿法(割线法)两大类：

牛顿法：原始牛顿法、修正牛顿法；
拟牛顿法：逆Broyden秩1法、逆Broyden秩1第二公式、BFS秩2法；

上面这5种方法的函数表达式、计算流程、代码都会详细说明和展示。并且，这5种方法都很经济(算的快)、实用(适用各种非刁钻的函数)、易用(计算流程好理解，便于编程)。

本文侧重的是方法的使用，不提推导和太具体的其他细节。

非线性方程组解法 —— 牛顿法

本文要介绍的5种方法可分为两大类：牛顿法、拟牛顿法。
先把两类方法的优缺点列出来：

牛顿法：用jacobi矩阵(导数)

优点：导数法收敛速度巨快(平方收敛)；自校正误差不会传递；
缺点：求导稍费事；只要赋值后的jacobi矩阵存在稀疏性、奇异性、病态等，就跪了；

拟牛顿法：割线法思想，用近似矩阵趋近jacobi矩阵

优点：jacobi矩阵的问题在这里都不是问题！这个优点极大提高解法的稳定性！！！
缺点：收敛速度介于平方收敛和直线收敛之间，稍慢一丢丢。

其实，牛顿法看上去要求导很麻烦，其实导数就求了一次而已！剩下都是在循环里对jacobi矩阵的赋值！所以去求导其实不是大问题。大问题就是每次赋值后的jacobi矩阵要求逆！！这就对数值矩阵的要求很高了！并且实际问题中经常出现赋值后的jacobi矩阵是稀疏的。
所以，如果原函数们不刁钻，两种方法都可胜任。但如果稍微感觉原函数有些复杂时，建议牺牲一丢的速度选择拟牛顿法！拟牛顿法的稳定性真的提高一个量级。

下面我们将正式开始方法的介绍，给定一个 $n\times n$ 的非线性方程组的通式：

$F(x) = [f_1(x),f_2(x),\cdots,f_n(x)]^{T} = 0$

对应的jacobi矩阵为：

$F^{'}(x) = \left( \begin{matrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_n} \\ \vdots & \vdots & & \vdots \\ \frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \cdots & \frac{\partial f_n}{\partial x_n} \end{matrix} \right)$

原始牛顿法

$x_0$ 自己给定，然后进行迭代：

$x^{k+1} = x^{k} - [F^{'}(x^{k})]^{-1}F(x^k)$

上述迭代可以用了，但是每次循环都要求逆很慢！所以换成下面这种形式：

$\begin{cases} x^{k+1} = x^{k} + \Delta x^{k} & 正常迭代\\ F^{'}(x^{k})\Delta x^{k} + F(x^k) = 0 & 每次解这个线性方程组 \end{cases}$

设真实解为 $x^{*}$ ，实现流程如下：

步1：给出初始近似 $x_0$ 及计算精度 $\varepsilon _1$ 和 $\varepsilon _2$ ；
步2：假定已进行了k次迭代，已求出 $x^{k}$ 和 $F(x^{k})$ ，计算 $F^{'}(x^{k})$ 并做赋值：
$A_k = F^{'}(x^{k}) \quad ；\quad b_k = F(x^{k})$
步3：解上面的线性方程组(用预处理后万能高斯-赛德尔迭代)：
$F^{'}(x^{k})\Delta x^{k} + F(x^k) = 0$
步4：求下面两个式子：
$x^{k+1} = x^{k} + \Delta x^{k} \quad ；\quad F(x^{k+1})$
步5：若 $\parallel \Delta x^{k} \parallel < \varepsilon _1$ 或 $\parallel F(x^{k+1}) \parallel < \varepsilon _2$ ，则置：
$x^{*} = x^{k+1}$
并转步6；否则进行如下赋值后回到步2：
$k = k+1 \quad ；\quad x^{k} = x^{k+1} \quad ；\quad F(x^{k}) = F(x^{k+1})$
步6：结束。

针对最开始的例子，下面给出matlab程序：

clear ; clc

% 未知数: 
syms x1 x2;

% 方程组: 统一用列向量
f1 = x1^2 - 10*x1 + x2^2 + 8;
f2 = x1*x2^2 + x1 - 10*x2 + 8;
% f1 = 4*x1 - x2 + 0.1*exp(x1)-1;
% f2 = -x1 + 4*x2 + 1/8*x1^2;
x = [x1;x2];
f = [f1;f2];

% 初始值: 统一用列向量
x0 = [2;3];
error_dxk = double( input('dxk范数的精度:') );
error_fkk = double( input('fkk范数的精度:') );
num = input('停止迭代次数:');

% jacobi1 = [diff(f1,x1) diff(f1,x2);diff(f2,x1) diff(f2,x2)]
% 直接用自带函数求雅克比矩阵:
jacobi = jacobian([f1,f2],[x1,x2]);

for k = 1:num
    Ak = double( subs(jacobi, x, x0) );
    bk = double( subs(f, x, x0) );
    dxk = pre_seidel(Ak,-bk,k);  % 步长
    x0 = x0 + dxk;
    fkk = double( subs(f, x, x0) );  % fk+1单纯用来判断
    if norm(dxk) < error_dxk | norm(fkk) < error_fkk
        break;
    end
end

if k < num
    x_result = x0;
    fprintf('精度已达要求，迭代提前结束!\n');
    fprintf('%d次迭代后, 近似解为:\n',k);
    x_result
else
    x_result = x0;
    fprintf('迭代次数已达上限!\n');
    fprintf('似解为:\n',k);
    x_result
end
    
fprintf('f1结果为:%f\n',double( subs(f1,x,x0) ));
fprintf('f2结果为:%f\n',double( subs(f2,x,x0) ));
fprintf('dxk范数:%f\n',norm(dxk));
fprintf('fkk范数:%f\n',norm(fkk));

结果：

dxk范数的精度:0.0001
fkk范数的精度:0.0001
停止迭代次数:200
第1次求解线性方程组, 当前万能赛德尔迭代矩阵谱半径为(越小越好): 0.8306
第2次求解线性方程组, 当前万能赛德尔迭代矩阵谱半径为(越小越好): 0.0668
第3次求解线性方程组, 当前万能赛德尔迭代矩阵谱半径为(越小越好): 0.1969
第4次求解线性方程组, 当前万能赛德尔迭代矩阵谱半径为(越小越好): 0.2209
第5次求解线性方程组, 当前万能赛德尔迭代矩阵谱半径为(越小越好): 0.2214
精度已达要求，迭代提前结束!
5次迭代后, 近似解为:

x_result =

    0.9999
    1.0000

f1结果为:0.000728
f2结果为:0.000184
dxk范数:0.000000
fkk范数:0.000751

说明：完全按照上面的流程编写的程序，非常好理解。
对应的pre_seidel自定义函数在下载地址里面都有。

修正牛顿法

对于上面的原始牛顿法，如果每步计算 $F^{'}(x^{k})$ 改为固定的 $F^{'}(x^{0})$ 可得：

$x^{k+1} = x^{k} - [F^{'}(x^{0})]^{-1}F(x^{k}) = x^{k} - BF(x^{k})$

这样就变成了"简化牛顿法"。很明显可以看到简化牛顿法是线性收敛的！
计算量小，但是收敛速度非常慢。

如果既拥有"原始牛顿法"的收敛速度，又拥有"简化牛顿法"的工作量省？—— 修正牛顿法。

对应的迭代格式为：

$\begin{cases} x^{k,0} = x^{k} & \\ x^{k,i} = x^{k,i-1} - [F^{'}(x^{k})]^{-1}\color{red}{F(x^{k,i-1})} & i=1,\cdots,m ; k = 1,\cdots,n\\ x^{k+1} = x^{k,m} & 条件3 \end{cases}$

从迭代公式可以得知，在每一个k的大循环内都有一个从1到m的小循环来求 $x^{k,m}$ ，下面同样给出实现流程：

步1：给出初始近似 $x_0$ 及计算精度 $\varepsilon _1$ 和 $\varepsilon _2$ ；
步2：根据方程阶数/个数n，求出效率最大化的m值或人为给定一个m值；
步3：假定已进行k此迭代，已算出 $x^{k}$ 和 $F(x^{k})$ ，计算并赋值：

$F^{'}(x^{k}) \quad ；\quad A_{k} = [F^{'}(x^{k})]^{-1}$

$\color{red}{步4(小循环)}$ ：进入每次的内层循环，假设已内层循环i次，已算出 $x^{k,i}$ 和 $F(x^{k,i})$ ，先做1个赋值，再做2个计算：

$b = F(x^{k,i}) \quad ；\quad x^{k,i+1} = x^{k,i} - A_kb \quad ；\quad F(x^{k,i+1})$

$\color{red}{步5(小循环)}$ ：先做3个赋值

$i = i+1 \quad ；\quad x^{k,i} = x^{k,i+1} \quad ；\quad F(x^{k,i}) = F(x^{k,i+1})$

然后做一个计算：

$z = A*b$

然后再一次判断：如果 $i == \color{red}{m}$ ，转到步6(出小循环)，否则回到步4(没出小循环)；

步6：若 $\parallel z \parallel < \varepsilon _1$ 或 $\parallel F(x^{k,\color{red}{m}}) \parallel < \varepsilon _2$ ，则赋值并结束：

$x^{*} = x^{k,m}$

否则， $k=k+1$ 然后重回步3(大循环)。

仍针对最开始的例子，给出matlab程序：

clear ; clc;

% 未知数: 
syms x1 x2;

% 方程组: 统一用列向量
f1 = x1^2 - 10*x1 + x2^2 + 8;
f2 = x1*x2^2 + x1 - 10*x2 + 8;
x = [x1;x2];
f = [f1;f2];

% 初始值: 统一用列向量
x0 = [5;4];
error_z = double( input('z范数的精度:') );
error_fk = double( input('fk范数的精度:') );
num = input('停止迭代次数:');

% 直接用自带函数求雅克比矩阵:
jacobi = jacobian([f1,f2],[x1,x2]);

% 小循环m取多少的判断: 和方程个数N有关,求w的最大值
syms M N;
mn0 = [N;M];
% 参数xzn是方程的个数:
xzn = length(x);  
% 原始的效率对比方程:
w = (N+1)*log(M+1)/( (N+M)*log(2) ); 
% 让w方程求最大值的xzm值:
xzm = double( solve( diff( subs(w,N,xzn),M ) ) );
mn1 = [xzn;xzm];
% 最高效率值:
wax = double( subs(w,mn0,mn1) );

% 开始修正牛顿迭代:
xki = x0;
for k = 1:num
    fk = double( subs(f,x,x0) );
    Ak = inv( double( subs(jacobi, x, x0) ) );
    for m = 1:round(xzm)
        b = fk;
        x0 = x0 - Ak*b;
        fki = double( subs(f,x,x0) );
        fk = fki;
        z = Ak*b;    
    end
    if norm(z) < error_z | norm(fk) < error_fk
        break;
    end
end

if k < num
    x_result = x0;
    fprintf('精度已达要求，迭代提前结束!\n');
    fprintf('%d次迭代后, 近似解为:\n',k);
    x_result
else
    x_result = x0;
    fprintf('迭代次数已达上限!\n');
    fprintf('似解为:\n',k);
    x_result
end

fprintf('f1结果为:%f\n',double( subs(f1,x,x0) ));
fprintf('f2结果为:%f\n',double( subs(f2,x,x0) ));
fprintf('z范数:%f\n',norm(z));
fprintf('fk范数:%f\n',norm(fk));

效果和原始牛顿法一样，可以感觉到速度提升了！！
这里用到了求最效率的m值，下面对它的求法加以补充(完全可以不求手动给)：

修正牛顿法 $N_{x}$ 与原始牛顿法 $N_{y}$ 的效率之比为：

$\frac{e(N_{x})}{e(N_{y})} = \frac{n+1}{n+m}\frac{ln^{(m+1)}}{ln^{2}}$

其中n就是方程组中方程的个数，这对于每个方程组都是固定常数。所以要想效率最高，那就让上式中右边含m参数的表达式取最大值即可(求导让导函数=0)，即如程序中所示。m求出若不是整数，就4舍5入。

至此，牛顿法的两种方法介绍完毕。下面简单对比一些二者的不同：

	原始牛顿法	修正牛顿法
区别	不求jacobi矩阵的逆每次要求线性方程组	不求线性方程组每次要求jacobi的逆

前文已经说明：每次要对矩阵求逆的话就很不稳定！
所以：其实原始牛顿法不错！好编程、速度快、不求逆、线性方程组有万能解法！修正牛顿法单纯提高了一丢丢效率，但是稳定性个人觉得变差了。
故，牛顿法中我推荐使用原始牛顿法。

非线性方程组解法 —— 拟牛顿法

牛顿法中要求导的jacobi矩阵，自然可以想到能否用差商来代替求导？
肯定是可以的，这就是割线法的思想，
即牛顿法是用超切平面趋近，割线法是用超割平面去趋近。
常用的割线法有：2点割线法、n点割线法、n+1点割线法；(n是方程个数)
其中n+1点割线法效率是最高的，但是稳定性没有n点割线法好！

所以，能不能构造一种算法，它具备n+1点割线法效率高的优点同时又增加稳定性？
这就是拟牛顿法。
所以：拟牛顿法是基于割线法，并做的比割线法更好的算法！
所以：拟牛顿法和割线法的"收敛速度"介于线性收敛和平方收敛之间；
所以：本文就不介绍割线法怎么搞，直接上最好的割线法 —— 拟牛顿法。

拟牛顿法中最好的是Broyden提出的操作，统称为Broyden方法。
思想：不断用一个低秩矩阵对 $A_k$ 进行修正；不同方法的区别就是那个低秩矩阵不同~
可以想象：低秩矩阵不同的秩数，就能得到一系列方法。
最常用：Broyden秩1、Broyden秩1第二方法、逆Broyden秩1，逆Broyden秩1第二方法。
另外还有秩2的方法，比如比较好的BFS。

秩1相关方法

1. Broyden秩1迭代公式：

$\begin{cases} x^{k+1} = x^{k} - A_{k}^{-1}F(x^{k}) & \\ A_{k+1} = A_{k} + \frac{ (y_{k} - A_{k}s_{k})(s_{k})^{T} }{ (s_{k})^{T}s_{k} } & \end{cases}$

其中： $s_{k} = x^{k+1} - x^{k} \quad ；\quad y_{k} = F(x^{k+1}) - F(x^{k})$

2. Broyden秩1第二方法迭代公式：

$\begin{cases} x^{k+1} = x^{k} - A_{k}^{-1}F(x^{k}) & \\ A_{k+1} = A_{k} + \frac{ F(x^{k+1}) F(x^{k+1})^{T} }{ F(x^{k+1})^{T} (x^{k+1} - x^{k}) } & \end{cases}$

上面的两种方法中的第一个式子都要求逆，不好看！我们用 $B_k = A_{k}^{-1}$ 代替，然后对上面的两个方法分别做改写，就可以得到它们的逆方法！

3. 逆Broyden秩1迭代公式：√

$\begin{cases} x^{k+1} = x^{k} - B_{k}F(x^{k}) & \\ B_{k+1} = B_{k} + \frac{ (s_{k} - B_{k}y_{k})(s_{k})^{T}B_{k} }{ (s_{k})^{T}B_{k}y_{k} } & \end{cases}$

4. 逆Broyden秩1第二方法迭代公式：√

$\begin{cases} x^{k+1} = x^{k} - B_{k}F(x^{k}) & \\ B_{k+1} = B_{k} + (s_{k} - B_{k}y_{k})\frac{ (s_{k} - B_{k}y_{k})^{T} }{(s_{k} - B_{k}y_{k})^{T}y_{k} } & \end{cases}$

对于上面4个方法/迭代公式，有下面3点说明：

文献中说"第二方法"类只适用于jacobi矩阵对称！但其实不对称的时候也可以用！不过稳定性大幅变差；
"逆方法"类的稳定性会提高！所以实际中/本文使用逆方法类。
方法间的区别只在第2个公式中右边那个项。故编程只用换那个表达式即可。

下面给出"逆Broyden秩1方法"实现流程：

步1：给出初始近似 $x_0$ 及计算精度 $\varepsilon _1$ 和 $\varepsilon _2$ ；
步2：计算初始矩阵 $B_{0}$ ，一般取：

$B_{0} = [F^{'}(x^{0})]^{-1}$

步3：k = 0；计算 $F(x^{0})$
步4：计算 $s_{k}$ 和 $x_{k+1}$ ：

$s_{k} = -B_{k}F(x^{k}) \quad ；\quad x^{k+1} = x^{k} + s_{k}$

步5：计算 $F(x^{k+1})$ ；检验若若 $\parallel s_{k} \parallel < \varepsilon _1$ 或 $\parallel F(x^{k+1}) \parallel < \varepsilon _2$ 则转步8，否则转步6；
步6：计算 $y_{k}$ ，并根据上面公式计算 $B_{k+1}$ ：

$y_{k} = F(x^{k+1}) - F(x^{k})$

$B_{k+1} = B_{k} + \frac{ (s_{k} - B_{k}y_{k})(s_{k})^{T}B_{k} }{ (s_{k})^{T}B_{k}y_{k} }$

步7：做4个赋值，然后回到步4：

$k = k+1 \quad ；\quad F(x^{k}) = F(x^{k+1}) \quad ；\quad B_{k} = B_{k+1} \quad ；\quad x^{k} = x^{k+1}$

步8： $x^{*} = x^{k+1}$ ，结束。

仍对于最开始的例子，给出逆Broyden秩1法的matlab程序：

clear; clc;

% 未知数: 
syms x1 x2;

% 方程组: 统一用列向量
f1 = x1^2 - 10*x1 + x2^2 + 8;
f2 = x1*x2^2 + x1 - 10*x2 + 8;
x = [x1;x2];
f = [f1;f2];

% 初始值: 统一用列向量
x0 = [2;3];
error_sk = double( input('sk范数的精度:') );
error_fkk = double( input('fkk范数的精度:') );
num = input('停止迭代次数:');

% 直接用自带函数求雅克比矩阵:
jacobi = jacobian([f1,f2],[x1,x2]);

% 开始求解前的初值:
Bk = inv( double( subs(jacobi,x,x0) ) );
fk = double( subs(f,x,x0) );
% 循环完全按照流程来
for k = 1:num
    sk = -Bk*fk;
    x0 = x0 + sk;
    fkk = double( subs(f,x,x0) );
    if norm(sk) < error_sk | norm(fkk) < error_fkk 
        break;
    end
    yk = fkk - fk;
    Bkk = Bk + (sk-Bk*yk)*(sk')*Bk/( sk'*Bk*yk );  % 不同方法改这里表达式即可
    fk = fkk;
    Bk = Bkk;
end

if k < num
    x_result = x0;
    fprintf('精度已达要求，迭代提前结束!\n');
    fprintf('%d次迭代后, 近似解为:\n',k);
    x_result
else
    x_result = x0;
    fprintf('迭代次数已达上限!\n');
    fprintf('似解为:\n',k);
    x_result
end

fprintf('f1结果为:%f\n',double( subs(f1,x,x0) ));
fprintf('f2结果为:%f\n',double( subs(f2,x,x0) ));
fprintf('sk范数:%f\n',norm(sk));
fprintf('fkk范数:%f\n',norm(fkk));

效果就是比牛顿法多迭代几次而已。
对应的逆Broyden秩1第二方法只用改程序中的Bkk表达式即可~ 故不多展示。

秩2相关方法

个人感觉：秩越高，"局部收敛"的范围更广一些(是好事)！后面会举例子说明这一点。

个人感觉稳定、效率高、效果好的秩2方法是：BFS(Broyden-Fletcher-Shanme)，迭代公式为：

$\begin{cases} x^{k+1} = x^{k} - B_{k}F(x^{k}) & \\ B_{k+1} = B_{k} + \frac{\mu _{k}s_{k}(s_{k})^{T} - s_{k}(y_{k})^{T}B_{k} - B_{k}y_{k}(s_{k})^{T}}{(s_{k})^{T}y_{k}}& \end{cases}$

其中 $\mu _{k}$ 的表达式为：

$\mu _{k} = 1 + \frac{(y_k)^{T}B_{k}y_{k}}{(s_{k})^{T}y_{k}}$

实现上和秩1法基本上就是一样的，只不过每次循环多算一个 $\mu _{k}$ 就行。给出程序：

clear; clc;

% 未知数: 
syms x1 x2;

% 方程组: 统一用列向量
f1 = x1^2 - 10*x1 + x2^2 + 8;
f2 = x1*x2^2 + x1 - 10*x2 + 8;
x = [x1;x2];
f = [f1;f2];

% 初始值: 统一用列向量
x0 = [5;4];
error_sk = double( input('sk范数的精度:') );
error_fkk = double( input('fkk范数的精度:') );
num = input('停止迭代次数:');

% 直接用自带函数求雅克比矩阵:
jacobi = jacobian([f1,f2],[x1,x2]);

% 开始求解前的初值:
Bk = inv( double( subs(jacobi,x,x0) ) );
fk = double( subs(f,x,x0) );

for k = 1:num
    sk = -Bk*fk;
    x0 = x0 + sk;
    fkk = double( subs(f,x,x0) );
    if norm(sk) < error_sk | norm(fkk) < error_fkk 
        break;
    end
    yk = fkk - fk;
    miuk = 1 + yk'*Bk*yk/(sk'*yk);   % 就多算一个这个而已
    Bkk = Bk + (miuk*sk*sk' - sk*yk'*Bk - Bk*yk*sk')/(sk'*yk);
    fk = fkk;
    Bk = Bkk;
end

if k < num
    x_result = x0;
    fprintf('精度已达要求，迭代提前结束!\n');
    fprintf('%d次迭代后, 近似解为:\n',k);
    x_result
else
    x_result = x0;
    fprintf('迭代次数已达上限!\n');
    fprintf('似解为:\n',k);
    x_result
end

fprintf('f1结果为:%f\n',double( subs(f1,x,x0) ));
fprintf('f2结果为:%f\n',double( subs(f2,x,x0) ));
fprintf('sk范数:%f\n',norm(sk));
fprintf('fkk范数:%f\n',norm(fkk));

至此，所有的方法就介绍完毕了。

总结

5种方法都在上面介绍并给出了程序，下面还有3点自己的心得体会(不一定正确！)。

1. 关于两套精度判断： $\varepsilon _1$ 和 $\varepsilon _2$ ：
不同方法衡量的对象稍有不同，大体可归纳为：
$\varepsilon _1$ 衡量是前后两次解的差值的范数； $\varepsilon _2$ 衡量是矩阵数值的范数；
可以认为是2套独立的衡量标准，只不过我们在线性方程组里常用的是第一种而已罢了；
只要迭代在收敛，这两个范数都是在不断缩小的往0趋近的！故其实都可设置成0.0001这种形式；
谁先到，意味着有一个衡量标准已经达到了，就可以结束了！
当然让两个都达到了再结束也可以。一句话：只要迭代在收敛，两套范数标准都是在下降趋向0的。

2. 拟牛顿法秩越高的方法"局部收敛"的范围会广一些！
对应上面同样的例子： $\varepsilon _1$ 和 $\varepsilon _2$ 都是0.0001；5种方法的初值都是[5;4]时：
只有BFS这个秩2方法的最终解是[1;1]
其他4种方法全找的是[2.1934;3.0205]
虽然这两个都是非线性方程组的解，但是为什么只有秩2的方法能搜到离起始点较远的解？
个人猜测是：秩2法的"局部收敛"范围比秩1法更广！
但不管是秩几，只要是牛顿法的大类(牛顿+割线+拟牛顿)，都是局部收敛的！

3. 使用推荐：
若方程组很大，看重求解速度：使用原始牛顿法；
若看重求解的稳定性：拟牛顿法BFS秩2；

补充说明：牛顿法最有可能出问题的地方就是jacobi矩阵每次赋完值之后的数值矩阵！稳定性太依赖这个导函数矩阵。拟牛顿法稳定性非常好！并且个人感觉：秩越高稳定性越好、收敛速度越趋于平方收敛！

本文所有程序下载地址

参考宝书

李庆扬, 莫孜中, 祁力群. 非线性方程组的数值解法[M]. 科学出版社, 2005.
范金燕, 袁亚湘. 非线性方程组数值方法[M]. 科学出版社, 2018.
李星. 数值分析[M]. 科学出版社, 2014.

数值分析：非线性方程组求解