前言

非線性方程組，顧名思義就是未知數(shù)的冪除了不是1陈莽，其他都有可能渤昌！線性方程組其實(shí)只是非常小的一類(lèi)，非線性方程組才是大類(lèi)走搁！也正因此非線性方程組包含各種各樣的方程形式独柑，所以它的解和對(duì)應(yīng)的求解方法不可能像線性方程組那樣完美，即都是局部收斂的朱盐。先給出一個(gè)直觀的非線性方程組例子：

$\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}$

個(gè)人對(duì)兩個(gè)問(wèn)題的理解：

1. 非線性方程組如果有解群嗤，一般都有很多解！如何理解：

把方程組的解看成是各個(gè)函數(shù)圖像的交點(diǎn)的兵琳。我們知道非線性方程組的各個(gè)函數(shù)就都是復(fù)雜曲線/面，甚至是高緯空間里的復(fù)雜東西骇径；線性方程組的各個(gè)函數(shù)就是最簡(jiǎn)單的直線/面躯肌！各個(gè)復(fù)雜函數(shù)圖像間的相交機(jī)會(huì)很多，并且只要相交破衔，就是多個(gè)交點(diǎn)(因?yàn)榻痪€清女、交面里有無(wú)數(shù)的交點(diǎn))，也就是有多個(gè)解~

可以想象晰筛，非線性方程組有多解是很平常的一件事~ 對(duì)于復(fù)雜的非線性函數(shù)沒(méi)解才不正常嫡丙！
可以想象拴袭，這些解是等價(jià)的！沒(méi)有說(shuō)是等級(jí)更高曙博，誰(shuí)等級(jí)低一些拥刻。都是解！因?yàn)椋褐灰墙飧赣荆?strong>只滿足一個(gè)條件：讓方程組中的各個(gè)方程=0般哼。所以無(wú)法用什么評(píng)判標(biāo)準(zhǔn)(比如范數(shù))來(lái)說(shuō)哪個(gè)解的等級(jí)高一些或者效果更好一些。
注意：這里的解等價(jià)和欠定線性方程組通解中的唯一極小范數(shù)解不一樣惠窄！可以想象二者的區(qū)別：非線性方程組中的解都是實(shí)打?qū)嵈嬖诘恼裘撸欢范ň€性方程組中除了特解，其他通解中的解說(shuō)存在也行杆融，說(shuō)不存在那就是因?yàn)榉匠虠l件(個(gè)數(shù))都不夠楞卡！這些是啥都行的通解和非線性方程組中實(shí)打?qū)嵈嬖诘慕饪隙ú荒鼙龋?/p>

這樣的話各個(gè)非線性方程組的局部收斂性就可以理解，即：空間中有很多解時(shí)脾歇，我每次只能找一個(gè)蒋腮，那我找誰(shuí)？找離我出發(fā)點(diǎn)最近的那個(gè)解唄介劫。所以不同的出發(fā)點(diǎn)徽惋，就有可能找到不同的解，這就是局部收斂性座韵。

2. 為什么不能用簡(jiǎn)單的替換先變成線性方程組求解险绘，然后再解每一個(gè)非線性方程？

意思是：每個(gè)方程中把所有和 $x_1$ 有關(guān)的用一個(gè)變量代替誉碴，所有 $x_2$ 有關(guān)的用一個(gè)變量代替宦棺，即方程1中用： $w_1 = x_{1}^{2} - 10*x_{1}$ ， $w_2 = x_{2}^{2} + 8$
但是很明顯方程2的第一項(xiàng) $x_{1}*x_{2}^{2}$ 兩個(gè)變量相乘黔帕，沒(méi)法用變量代替~
并且代咸，即使在方程2中能代替，那么就會(huì)有 $w_3$ 和 $w_4$ 成黄，這樣總未知數(shù)變成4個(gè)而方程只有2個(gè)呐芥，還是解不了。
所以奋岁，非線性方程組不可能用簡(jiǎn)單的線性變量代換來(lái)解聂沙。

本文介紹最常用鹉戚、最好用的非線性方程組解法，包括牛頓法和擬牛頓法(割線法)兩大類(lèi)：

牛頓法：原始牛頓法、修正牛頓法缸匪；
擬牛頓法：逆Broyden秩1法、逆Broyden秩1第二公式、BFS秩2法；

上面這5種方法的函數(shù)表達(dá)式女嘲、計(jì)算流程、代碼都會(huì)詳細(xì)說(shuō)明和展示诞帐。并且欣尼，這5種方法都很經(jīng)濟(jì)(算的快)、實(shí)用(適用各種非刁鉆的函數(shù))景埃、易用(計(jì)算流程好理解媒至，便于編程)。

本文側(cè)重的是方法的使用谷徙，不提推導(dǎo)和太具體的其他細(xì)節(jié)拒啰。

非線性方程組解法 —— 牛頓法

本文要介紹的5種方法可分為兩大類(lèi)：牛頓法、擬牛頓法完慧。
先把兩類(lèi)方法的優(yōu)缺點(diǎn)列出來(lái)：

牛頓法：用jacobi矩陣(導(dǎo)數(shù))

優(yōu)點(diǎn)：導(dǎo)數(shù)法收斂速度巨快(平方收斂)谋旦；自校正誤差不會(huì)傳遞；
缺點(diǎn)：求導(dǎo)稍費(fèi)事屈尼；只要賦值后的jacobi矩陣存在稀疏性册着、奇異性、病態(tài)等脾歧，就跪了甲捏；

擬牛頓法：割線法思想，用近似矩陣趨近jacobi矩陣

優(yōu)點(diǎn)：jacobi矩陣的問(wèn)題在這里都不是問(wèn)題鞭执！這個(gè)優(yōu)點(diǎn)極大提高解法的穩(wěn)定性Ｋ径佟！兄纺！
缺點(diǎn)：收斂速度介于平方收斂和直線收斂之間大溜，稍慢一丟丟。

其實(shí)估脆，牛頓法看上去要求導(dǎo)很麻煩钦奋，其實(shí)導(dǎo)數(shù)就求了一次而已！剩下都是在循環(huán)里對(duì)jacobi矩陣的賦值疙赠！所以去求導(dǎo)其實(shí)不是大問(wèn)題付材。大問(wèn)題就是每次賦值后的jacobi矩陣要求逆！圃阳！這就對(duì)數(shù)值矩陣的要求很高了伞租！并且實(shí)際問(wèn)題中經(jīng)常出現(xiàn)賦值后的jacobi矩陣是稀疏的。
所以限佩，如果原函數(shù)們不刁鉆，兩種方法都可勝任。但如果稍微感覺(jué)原函數(shù)有些復(fù)雜時(shí)祟同，建議犧牲一丟的速度選擇擬牛頓法作喘！擬牛頓法的穩(wěn)定性真的提高一個(gè)量級(jí)。

下面我們將正式開(kāi)始方法的介紹晕城，給定一個(gè) $n\times n$ 的非線性方程組的通式：

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

對(duì)應(yīng)的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$ 自己給定泞坦，然后進(jìn)行迭代：

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

上述迭代可以用了，但是每次循環(huán)都要求逆很慢砖顷！所以換成下面這種形式：

$\begin{cases} x^{k+1} = x^{k} + \Delta x^{k} & 正常迭代\\ F^{'}(x^{k})\Delta x^{k} + F(x^k) = 0 & 每次解這個(gè)線性方程組 \end{cases}$

設(shè)真實(shí)解為 $x^{*}$ 贰锁，實(shí)現(xiàn)流程如下：

步1：給出初始近似 $x_0$ 及計(jì)算精度 $\varepsilon _1$ 和 $\varepsilon _2$ ；
步2：假定已進(jìn)行了k次迭代滤蝠，已求出 $x^{k}$ 和 $F(x^{k})$ 豌熄，計(jì)算 $F^{'}(x^{k})$ 并做賦值：
$A_k = F^{'}(x^{k}) \quad ；\quad b_k = F(x^{k})$
步3：解上面的線性方程組(用預(yù)處理后萬(wàn)能高斯-賽德?tīng)柕?/a>)：
$F^{'}(x^{k})\Delta x^{k} + F(x^k) = 0$

步4：求下面兩個(gè)式子：
$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}$
并轉(zhuǎn)步6；否則進(jìn)行如下賦值后回到步2：
$k = k+1 \quad 览闰；\quad x^{k} = x^{k+1} \quad 芯肤；\quad F(x^{k}) = F(x^{k+1})$

步6：結(jié)束。

針對(duì)最開(kāi)始的例子压鉴，下面給出matlab程序：

clear ; clc

% 未知數(shù): 
syms x1 x2;

% 方程組: 統(tǒng)一用列向量
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];

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

% jacobi1 = [diff(f1,x1) diff(f1,x2);diff(f2,x1) diff(f2,x2)]
% 直接用自帶函數(shù)求雅克比矩陣:
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);  % 步長(zhǎng)
    x0 = x0 + dxk;
    fkk = double( subs(f, x, x0) );  % fk+1單純用來(lái)判斷
    if norm(dxk) < error_dxk | norm(fkk) < error_fkk
        break;
    end
end

if k < num
    x_result = x0;
    fprintf('精度已達(dá)要求崖咨，迭代提前結(jié)束!\n');
    fprintf('%d次迭代后, 近似解為:\n',k);
    x_result
else
    x_result = x0;
    fprintf('迭代次數(shù)已達(dá)上限!\n');
    fprintf('似解為:\n',k);
    x_result
end
    
fprintf('f1結(jié)果為:%f\n',double( subs(f1,x,x0) ));
fprintf('f2結(jié)果為:%f\n',double( subs(f2,x,x0) ));
fprintf('dxk范數(shù):%f\n',norm(dxk));
fprintf('fkk范數(shù):%f\n',norm(fkk));

結(jié)果：

dxk范數(shù)的精度:0.0001
fkk范數(shù)的精度:0.0001
停止迭代次數(shù):200
第1次求解線性方程組, 當(dāng)前萬(wàn)能賽德?tīng)柕仃囎V半徑為(越小越好): 0.8306
第2次求解線性方程組, 當(dāng)前萬(wàn)能賽德?tīng)柕仃囎V半徑為(越小越好): 0.0668
第3次求解線性方程組, 當(dāng)前萬(wàn)能賽德?tīng)柕仃囎V半徑為(越小越好): 0.1969
第4次求解線性方程組, 當(dāng)前萬(wàn)能賽德?tīng)柕仃囎V半徑為(越小越好): 0.2209
第5次求解線性方程組, 當(dāng)前萬(wàn)能賽德?tīng)柕仃囎V半徑為(越小越好): 0.2214
精度已達(dá)要求，迭代提前結(jié)束!
5次迭代后, 近似解為:

x_result =

    0.9999
    1.0000

f1結(jié)果為:0.000728
f2結(jié)果為:0.000184
dxk范數(shù):0.000000
fkk范數(shù):0.000751

說(shuō)明：完全按照上面的流程編寫(xiě)的程序油吭，非常好理解击蹲。
對(duì)應(yīng)的pre_seidel自定義函數(shù)在下載地址里面都有。

修正牛頓法

對(duì)于上面的原始牛頓法上鞠，如果每步計(jì)算 $F^{'}(x^{k})$ 改為固定的 $F^{'}(x^{0})$ 可得：

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

這樣就變成了"簡(jiǎn)化牛頓法"际邻。很明顯可以看到簡(jiǎn)化牛頓法是線性收斂的！
計(jì)算量小芍阎，但是收斂速度非常慢世曾。

如果既擁有"原始牛頓法"的收斂速度，又擁有"簡(jiǎn)化牛頓法"的工作量是聪獭轮听？—— 修正牛頓法。

對(duì)應(yīng)的迭代格式為：

$\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}$

從迭代公式可以得知岭佳，在每一個(gè)k的大循環(huán)內(nèi)都有一個(gè)從1到m的小循環(huán)來(lái)求 $x^{k,m}$ 血巍，下面同樣給出實(shí)現(xiàn)流程：

步1：給出初始近似 $x_0$ 及計(jì)算精度 $\varepsilon _1$ 和 $\varepsilon _2$ ；
步2：根據(jù)方程階數(shù)/個(gè)數(shù)n珊随，求出效率最大化的m值或人為給定一個(gè)m值述寡；
步3：假定已進(jìn)行k此迭代柿隙，已算出 $x^{k}$ 和 $F(x^{k})$ ，計(jì)算并賦值：

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

$\color{red}{步4(小循環(huán))}$ ：進(jìn)入每次的內(nèi)層循環(huán)禀崖，假設(shè)已內(nèi)層循環(huán)i次，已算出 $x^{k,i}$ 和 $F(x^{k,i})$ 螟炫，先做1個(gè)賦值波附，再做2個(gè)計(jì)算：

$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(小循環(huán))}$ ：先做3個(gè)賦值

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

然后做一個(gè)計(jì)算：

$z = A*b$

然后再一次判斷：如果 $i == \color{red}{m}$ 然评，轉(zhuǎn)到步6(出小循環(huán))仅财，否則回到步4(沒(méi)出小循環(huán))；

步6：若 $\parallel z \parallel < \varepsilon _1$ 或 $\parallel F(x^{k,\color{red}{m}}) \parallel < \varepsilon _2$ 沾瓦，則賦值并結(jié)束：

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

否則满着， $k=k+1$ 然后重回步3(大循環(huán))。

仍針對(duì)最開(kāi)始的例子贯莺，給出matlab程序：

clear ; clc;

% 未知數(shù): 
syms x1 x2;

% 方程組: 統(tǒng)一用列向量
f1 = x1^2 - 10*x1 + x2^2 + 8;
f2 = x1*x2^2 + x1 - 10*x2 + 8;
x = [x1;x2];
f = [f1;f2];

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

% 直接用自帶函數(shù)求雅克比矩陣:
jacobi = jacobian([f1,f2],[x1,x2]);

% 小循環(huán)m取多少的判斷: 和方程個(gè)數(shù)N有關(guān),求w的最大值
syms M N;
mn0 = [N;M];
% 參數(shù)xzn是方程的個(gè)數(shù):
xzn = length(x);  
% 原始的效率對(duì)比方程:
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) );

% 開(kāi)始修正牛頓迭代:
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('精度已達(dá)要求风喇，迭代提前結(jié)束!\n');
    fprintf('%d次迭代后, 近似解為:\n',k);
    x_result
else
    x_result = x0;
    fprintf('迭代次數(shù)已達(dá)上限!\n');
    fprintf('似解為:\n',k);
    x_result
end

fprintf('f1結(jié)果為:%f\n',double( subs(f1,x,x0) ));
fprintf('f2結(jié)果為:%f\n',double( subs(f2,x,x0) ));
fprintf('z范數(shù):%f\n',norm(z));
fprintf('fk范數(shù):%f\n',norm(fk));

效果和原始牛頓法一樣，可以感覺(jué)到速度提升了Ｂ铺健魂莫！
這里用到了求最效率的m值，下面對(duì)它的求法加以補(bǔ)充(完全可以不求手動(dòng)給)：

修正牛頓法 $N_{x}$ 與原始牛頓法 $N_{y}$ 的效率之比為：

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

其中n就是方程組中方程的個(gè)數(shù)爹耗，這對(duì)于每個(gè)方程組都是固定常數(shù)耙考。所以要想效率最高，那就讓上式中右邊含m參數(shù)的表達(dá)式取最大值即可(求導(dǎo)讓導(dǎo)函數(shù)=0)潭兽，即如程序中所示倦始。m求出若不是整數(shù)，就4舍5入山卦。

至此鞋邑，牛頓法的兩種方法介紹完畢。下面簡(jiǎn)單對(duì)比一些二者的不同：

?	原始牛頓法	修正牛頓法
區(qū)別	不求jacobi矩陣的逆每次要求線性方程組	不求線性方程組每次要求jacobi的逆

前文已經(jīng)說(shuō)明：每次要對(duì)矩陣求逆的話就很不穩(wěn)定账蓉！
所以：其實(shí)原始牛頓法不錯(cuò)枚碗！好編程、速度快铸本、不求逆肮雨、線性方程組有萬(wàn)能解法！修正牛頓法單純提高了一丟丟效率箱玷，但是穩(wěn)定性個(gè)人覺(jué)得變差了怨规。
故陌宿，牛頓法中我推薦使用原始牛頓法。

非線性方程組解法 —— 擬牛頓法

牛頓法中要求導(dǎo)的jacobi矩陣椅亚，自然可以想到能否用差商來(lái)代替求導(dǎo)限番？
肯定是可以的，這就是割線法的思想呀舔，
即牛頓法是用超切平面趨近，割線法是用超割平面去趨近扩灯。
常用的割線法有：2點(diǎn)割線法媚赖、n點(diǎn)割線法、n+1點(diǎn)割線法珠插；(n是方程個(gè)數(shù))
其中n+1點(diǎn)割線法效率是最高的惧磺，但是穩(wěn)定性沒(méi)有n點(diǎn)割線法好！

所以捻撑，能不能構(gòu)造一種算法磨隘，它具備n+1點(diǎn)割線法效率高的優(yōu)點(diǎn)同時(shí)又增加穩(wěn)定性？
這就是擬牛頓法顾患。
所以：擬牛頓法是基于割線法番捂，并做的比割線法更好的算法！
所以：擬牛頓法和割線法的"收斂速度"介于線性收斂和平方收斂之間江解；
所以：本文就不介紹割線法怎么搞设预，直接上最好的割線法 —— 擬牛頓法。

擬牛頓法中最好的是Broyden提出的操作犁河，統(tǒng)稱(chēng)為Broyden方法鳖枕。
思想：不斷用一個(gè)低秩矩陣對(duì) $A_k$ 進(jìn)行修正；不同方法的區(qū)別就是那個(gè)低秩矩陣不同~
可以想象：低秩矩陣不同的秩數(shù)桨螺，就能得到一系列方法宾符。
最常用：Broyden秩1、Broyden秩1第二方法灭翔、逆Broyden秩1魏烫，逆Broyden秩1第二方法。
另外還有秩2的方法缠局，比如比較好的BFS则奥。

秩1相關(guān)方法

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}$

上面的兩種方法中的第一個(gè)式子都要求逆狭园，不好看读处！我們用 $B_k = A_{k}^{-1}$ 代替，然后對(duì)上面的兩個(gè)方法分別做改寫(xiě)唱矛，就可以得到它們的逆方法罚舱！

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}$

對(duì)于上面4個(gè)方法/迭代公式井辜，有下面3點(diǎn)說(shuō)明：

文獻(xiàn)中說(shuō)"第二方法"類(lèi)只適用于jacobi矩陣對(duì)稱(chēng)！但其實(shí)不對(duì)稱(chēng)的時(shí)候也可以用管闷！不過(guò)穩(wěn)定性大幅變差粥脚；
"逆方法"類(lèi)的穩(wěn)定性會(huì)提高！所以實(shí)際中/本文使用逆方法類(lèi)包个。
方法間的區(qū)別只在第2個(gè)公式中右邊那個(gè)項(xiàng)刷允。故編程只用換那個(gè)表達(dá)式即可。

下面給出"逆Broyden秩1方法"實(shí)現(xiàn)流程：

步1：給出初始近似 $x_0$ 及計(jì)算精度 $\varepsilon _1$ 和 $\varepsilon _2$ 碧囊；
步2：計(jì)算初始矩陣 $B_{0}$ 树灶，一般取：

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

步3：k = 0糯而；計(jì)算 $F(x^{0})$
步4：計(jì)算 $s_{k}$ 和 $x_{k+1}$ ：

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

步5：計(jì)算 $F(x^{k+1})$ ；檢驗(yàn)若若 $\parallel s_{k} \parallel < \varepsilon _1$ 或 $\parallel F(x^{k+1}) \parallel < \varepsilon _2$ 則轉(zhuǎn)步8熄驼，否則轉(zhuǎn)步6像寒；
步6：計(jì)算 $y_{k}$ ，并根據(jù)上面公式計(jì)算 $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個(gè)賦值瓜贾，然后回到步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}$ 序臂，結(jié)束。

仍對(duì)于最開(kāi)始的例子实束，給出逆Broyden秩1法的matlab程序：

clear; clc;

% 未知數(shù): 
syms x1 x2;

% 方程組: 統(tǒng)一用列向量
f1 = x1^2 - 10*x1 + x2^2 + 8;
f2 = x1*x2^2 + x1 - 10*x2 + 8;
x = [x1;x2];
f = [f1;f2];

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

% 直接用自帶函數(shù)求雅克比矩陣:
jacobi = jacobian([f1,f2],[x1,x2]);

% 開(kāi)始求解前的初值:
Bk = inv( double( subs(jacobi,x,x0) ) );
fk = double( subs(f,x,x0) );
% 循環(huán)完全按照流程來(lái)
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 );  % 不同方法改這里表達(dá)式即可
    fk = fkk;
    Bk = Bkk;
end

if k < num
    x_result = x0;
    fprintf('精度已達(dá)要求奥秆，迭代提前結(jié)束!\n');
    fprintf('%d次迭代后, 近似解為:\n',k);
    x_result
else
    x_result = x0;
    fprintf('迭代次數(shù)已達(dá)上限!\n');
    fprintf('似解為:\n',k);
    x_result
end

fprintf('f1結(jié)果為:%f\n',double( subs(f1,x,x0) ));
fprintf('f2結(jié)果為:%f\n',double( subs(f2,x,x0) ));
fprintf('sk范數(shù):%f\n',norm(sk));
fprintf('fkk范數(shù):%f\n',norm(fkk));

效果就是比牛頓法多迭代幾次而已。
對(duì)應(yīng)的逆Broyden秩1第二方法只用改程序中的Bkk表達(dá)式即可~ 故不多展示咸灿。

秩2相關(guān)方法

個(gè)人感覺(jué)：秩越高构订，"局部收斂"的范圍更廣一些(是好事)！后面會(huì)舉例子說(shuō)明這一點(diǎn)避矢。

個(gè)人感覺(jué)穩(wěn)定悼瘾、效率高、效果好的秩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}$ 的表達(dá)式為：

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

實(shí)現(xiàn)上和秩1法基本上就是一樣的亥宿，只不過(guò)每次循環(huán)多算一個(gè) $\mu _{k}$ 就行。給出程序：

clear; clc;

% 未知數(shù): 
syms x1 x2;

% 方程組: 統(tǒng)一用列向量
f1 = x1^2 - 10*x1 + x2^2 + 8;
f2 = x1*x2^2 + x1 - 10*x2 + 8;
x = [x1;x2];
f = [f1;f2];

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

% 直接用自帶函數(shù)求雅克比矩陣:
jacobi = jacobian([f1,f2],[x1,x2]);

% 開(kāi)始求解前的初值:
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);   % 就多算一個(gè)這個(gè)而已
    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('精度已達(dá)要求砂沛，迭代提前結(jié)束!\n');
    fprintf('%d次迭代后, 近似解為:\n',k);
    x_result
else
    x_result = x0;
    fprintf('迭代次數(shù)已達(dá)上限!\n');
    fprintf('似解為:\n',k);
    x_result
end

fprintf('f1結(jié)果為:%f\n',double( subs(f1,x,x0) ));
fprintf('f2結(jié)果為:%f\n',double( subs(f2,x,x0) ));
fprintf('sk范數(shù):%f\n',norm(sk));
fprintf('fkk范數(shù):%f\n',norm(fkk));

至此烫扼，所有的方法就介紹完畢了。

總結(jié)

5種方法都在上面介紹并給出了程序碍庵，下面還有3點(diǎn)自己的心得體會(huì)(不一定正確映企！)悟狱。

1. 關(guān)于兩套精度判斷： $\varepsilon _1$ 和 $\varepsilon _2$ ：
不同方法衡量的對(duì)象稍有不同，大體可歸納為：
$\varepsilon _1$ 衡量是前后兩次解的差值的范數(shù)堰氓； $\varepsilon _2$ 衡量是矩陣數(shù)值的范數(shù)挤渐；
可以認(rèn)為是2套獨(dú)立的衡量標(biāo)準(zhǔn)，只不過(guò)我們?cè)诰€性方程組里常用的是第一種而已罷了双絮；
只要迭代在收斂浴麻，這兩個(gè)范數(shù)都是在不斷縮小的往0趨近的！故其實(shí)都可設(shè)置成0.0001這種形式掷邦；
誰(shuí)先到白胀，意味著有一個(gè)衡量標(biāo)準(zhǔn)已經(jīng)達(dá)到了，就可以結(jié)束了抚岗！
當(dāng)然讓兩個(gè)都達(dá)到了再結(jié)束也可以。一句話：只要迭代在收斂哪怔，兩套范數(shù)標(biāo)準(zhǔn)都是在下降趨向0的宣蔚。

2. 擬牛頓法秩越高的方法"局部收斂"的范圍會(huì)廣一些！
對(duì)應(yīng)上面同樣的例子： $\varepsilon _1$ 和 $\varepsilon _2$ 都是0.0001认境；5種方法的初值都是[5;4]時(shí)：
只有BFS這個(gè)秩2方法的最終解是[1;1]
其他4種方法全找的是[2.1934;3.0205]
雖然這兩個(gè)都是非線性方程組的解胚委，但是為什么只有秩2的方法能搜到離起始點(diǎn)較遠(yuǎn)的解？
個(gè)人猜測(cè)是：秩2法的"局部收斂"范圍比秩1法更廣叉信！
但不管是秩幾亩冬，只要是牛頓法的大類(lèi)(牛頓+割線+擬牛頓)，都是局部收斂的硼身！

3. 使用推薦：
若方程組很大硅急，看重求解速度：使用原始牛頓法；
若看重求解的穩(wěn)定性：擬牛頓法BFS秩2佳遂；

補(bǔ)充說(shuō)明：牛頓法最有可能出問(wèn)題的地方就是jacobi矩陣每次賦完值之后的數(shù)值矩陣营袜！穩(wěn)定性太依賴(lài)這個(gè)導(dǎo)函數(shù)矩陣。擬牛頓法穩(wěn)定性非常好丑罪！并且個(gè)人感覺(jué)：秩越高穩(wěn)定性越好荚板、收斂速度越趨于平方收斂！

本文所有程序下載地址

參考寶書(shū)

李慶揚(yáng), 莫孜中, 祁力群. 非線性方程組的數(shù)值解法[M]. 科學(xué)出版社, 2005.
范金燕, 袁亞湘. 非線性方程組數(shù)值方法[M]. 科學(xué)出版社, 2018.
李星. 數(shù)值分析[M]. 科學(xué)出版社, 2014.

數(shù)值分析：非線性方程組求解