在matlab中使用牛頓法求解非線性方程組

引言

我之前因為需要求解非線性方程組，在各種參考書上尋找現成的方程組求解算例，但發現大部分的書中的程序都是有問題的，要不然是輸出出來的結果精度不夠，要不然就根本跑不了。所以我自己參考了數值分析方法[1],非線性數值分析的理論與方法[2],還有一些單自變數牛頓法求解方程的算例，寫了下面的程序，如果不需要了解牛頓法的意義，可以直接從標題一開始看。

對於很多非線性方程組來說，數值解法是求解方程的唯一方法。在各類參考書中，對於牛頓法的編程介紹通常都是在一元方程中進行求解，很少有關於多自變數，非線性方程組的求解對於非線性方程組來說，方程的形式一般如下所示：

$m{y_i}=f(m{x_i}) space (i=1,2,...)=0$

進行數值迭代時，先選取一個初始點：

$m{x_0}=left(egin{array}{ccc} x_{01}\ x_{02}\x_{03}\ x_{04}\ vdots end{array} ight)$

將這個初始點代入已知的方程組中，便可得到初始點對應的函數值，假設函數的下一個迭代點滿足以下的關係：

當非奇，可以得到 $m{x^{k+1}}=m{x^k}-m{F(x^k)^{-1}}m{F(x^k)}(k=0,1,cdots)$ 其中 $m{F^{prime}left(x^{k} ight)}=left{egin{array}{cccc}{frac{partial f_{1}left(x^{k} ight)}{partial x_{1}}} & {frac{partial f_{1}left(x^{k} ight)}{partial x_{2}}} & {cdots} & {frac{partial f_{1}left(x^{k} ight)}{partial x_{n}}} \ {frac{partial f_{2}left(x^{k} ight)}{partial x_{1}}} & {frac{partial f_{2}left(x^{k} ight)}{partial x_{2}}} & {dots} & {frac{partial f_{2}left(x^{k} ight)}{partial x_{n}}} \ {ldots} & {ldots} & {ldots} & {ldots} \ {frac{partial f_{n}left(x^{k} ight)}{partial x_{1}}} & {frac{partial f_{n}left(x^{k} ight)}{partial x_{2}}} & {cdots} & {frac{partial f_{n}left(x^{k} ight)}{partial x_{n}}}end{array} ight}$

是在處的雅可比矩陣，通過初始點計算第二個迭代點，再通過第二個迭代點，計算第三個，以此類推，知道計算出滿足精度要求的點兩個迭代點之間的差值可以由如下公式計算得出：

$F^{prime}left(x^{k} ight) Delta x^{k}=-Fleft(x^{k} ight)$

這個方程就是牛頓方程，求得解 ,即可得到：

$x^{k+1}=x^{k}+Delta x^{k}$

再繼續迭代即可，當滿足精度要求，即 $|x^{k+1}-x^k|<epsilon$ 時停止迭代即可。

具體演算法在matlab裏可以這樣實現：

1編寫計算方程組值的函數

我所用的 matlab 版本為 R2015b。

首先假設要解的非線性方程組是如下形式：

$left{egin{array}{cccc} f_1(x_1,x_2)&=0.5sin x_1+0.1cos x_1x_2-x_1=0\ f_2(x_1,x_2)&=0.5cos x_1-0.1cos x_2-x_2=0 end{array} ight.$

然後建立一個函數文件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是滿足精度要求的最終迭代點。