引言
我之前因為需要求解非線性方程組，在各種參考書上尋找現成的方程組求解算例，但發現大部分的書中的程式都是有問題的，要不然是輸出出來的結果精度不夠，要不然就根本跑不了。是以我自己參考了
數值分析方法[1],
非線性數值分析的理論與方法[2],還有一些單自變量牛頓法求解方程的算例，寫了下面的程式，如果不需要了解牛頓法的意義，可以直接從
标題一開始看。
對于很多
非線性方程組來說，數值解法是求解方程的唯一方法。在各類參考書中，對于牛頓法的程式設計介紹通常都是在一進制方程中進行求解，很少有關于多自變量，
非線性方程組的求解 對于非線性方程組來說，方程的形式一般如下所示：
進行數值疊代時，先選取一個初始點 ：
将這個初始點代入已知的方程組中，便可得到初始點對應的函數值，假設函數的下一個疊代點滿足以下的關系：
當
非奇，可以得到
其中
是在
處的雅可比矩陣，通過初始點計算第二個疊代點，再通過第二個疊代點，計算第三個，以此類推，知道計算出滿足精度要求的點 兩個疊代點之間的內插補點可以由如下公式計算得出：
這個方程就是牛頓方程，求得解
,即可得到：
再繼續疊代即可，當滿足精度要求，即
時停止疊代即可。
具體算法在matlab裡可以這樣實作：
1編寫計算方程組值的函數
我所用的 matlab 版本為 R2015b。
首先假設要解的非線性方程組是如下形式 ：
然後建立一個函數檔案
fun.m
用于存放需要求解的非線性方程組。
function f = fun(x)
k=2;
for i=1:k
x(i)=sym (['x',num2str(i)]);
end
f(1)=0.5*sin(x(1))+0.1*cos(x(1)*x(2))-x(1);
f(2)=0.5*cos(x(1))-0.1*cos(x(2))-x(2);
end
在這裡k是所求解方程的個數，同時也是自變量的個數，第一個
for
循環是将自變量
,
變成符号，便于在數值計算函數中進行雅可比矩陣的計算，當輸入一個二進制數組時即可計算出來
和
。
2編寫算法
按照之前的牛頓法的分析，建立
mulNewton.m
的函數檔案：
function [allx,ally,r,n]=mulNewton(F,x0,eps)
if nargin==2
eps=1.0e-4;
end
x0 = transpose(x0);
Fx = subs(F,transpose(symvar(F)),x0);
var = transpose(symvar(F));
dF = jacobian(F,var);
dFx = subs(dF,transpose(symvar(F)),x0);
n=dFx;
r=x0-inv(dFx)*Fx';
n=1;
tol=1;
N=100;
symx=length(x0);
ally=zeros(symx,N);
allx=zeros(symx,N);
while tol>eps
x0=r;
Fx = subs(F,transpose(symvar(F)),x0);
dFx = subs(dF,transpose(symvar(F)),x0);
r=vpa(x0-inv(dFx)*Fx');
tol=norm(r-x0)
if(n>N)
disp('疊代步數太多，可能不收斂！');
break;
end
allx(:,n)=x0;
ally(:,n)=Fx;
n=n+1;
end
輸入的自變量
是需要求解的方程組檔案，也就是我們之前定義的
fun.m
檔案，
x_0
是一個數組，也就是我們選取的初始點，
eps
是可以不輸入的精度值，如果不輸入，就預設為1.0e-4 輸出值當中，
allx
是用來記錄每一步疊代的點的矩陣，
ally
是用來記錄疊代點對應計算出來的函數值的，
r
是滿足精度要求的最終疊代點。
3計算算例
首先先将
fun.m
和
mulNewton.m
檔案放到同一檔案夾下 然後，選取初始疊代點(1,2),在指令行中輸入如下指令：
即初始點選擇為
,精度為預設精度
,得到最終疊代結果如下所示：
即疊代4次，算出滿足精度要求的點為 (0.198,0.398)。
再看看各步疊代得到的
和
：
疊代到最後一步時，
值為(-1.56e-5,-3.545e-5),兩個值都十分接近零，說明計算得到的
點的精度還是比較令人滿意的。
4總結
對于希望解維數更多的非線性方程的同學來說，隻需要将
fun.m
檔案中的
更改為所需維數，檔案中的方程替換為你自己所需的方程即可。
參考文獻
[1]王能超. 數值分析簡明教程[M]. 2012.
[2]黃象鼎. 非線性數值分析的理論與方法[M]. 武漢大學出版社, 2004.