在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是满足精度要求的最终迭代点。