[转载]非线性方程组的数值解法
2013-03-08 19:03阅读：
非线性方程组的数值解法
摘要：本文介绍了几种求解非线性方程组的数值解法，并讨论了其对应的收敛速度.而且将其进行比较，阐述了各种数值解法的优缺点.而且列举了其对应的例子，将其进行解决.
0 、引言
n个变量n个方程（n>1）的方程组表示为：
（1）
为了求解非线性方程组的值，为此探讨几种求解非线性方程组的数值解法.而且列举了其对应的例子，将其进行解决.为此探讨如下：
一、基本概念及解法
1 、非线性方程组的一般解法
如果用矢量表示法，也可以表示为：
f(x)=0
(2)
在数值迭代中的一种等价的表示方法：x(k+1)=G(x(k))，k=0,1,……
(3)
由(2)和(3)的关系，就可以发现：f(x)=x-G(x)=0
(4)
数值迭代解法已有很多，但其实质不外乎是选定一个合适的x(k)作为起始点的解，然后用迭代方法导出一个修正的近似解x(k+1)使（1）式更加接近0.如果代入（1）式后方程得到满足则结束.如果（1）式还得不到满足，则必须再迭代一次，以便找到更加靠近答案的一组x(k+2).
这种方法算法为：
（1）
选择一个解的初始值x(0)；
（2）
用（3）式迭代取得一个较好的近似解；
（3）
当过程满足终止判定时停止迭代.
2 、直接迭代法
迭代法中最基本的方法为图形法和二分法.先从非线性方程开始，图形法根据已知方程用matlab画出其图形，从曲线与轴的交点可以看出其有多少个交点，就是对应方程有多少个解.二分法就是确定区间（a,b）后，取（a,b）的中点x(0)=(a+b)/2，若f(x(0))=0，则x(0)是根，否则，如f(a)*f(x(0))<0，令a1=a,b1=x(0)；如f(x(0))*f(b)<0，令a1=x(0),b1=b在(a1,b1)内至少有一个根，再取的中点如此进行下去，足够大时即可达到满意的精度.并由此推广到求解非线性方程组中去.
不过采用迭代法要考虑其收敛性，对于一维迭代法收敛的充分条件（非必要条件）是：
G(x)’=|dG(x)/dx|<1
如果此条件不满足，则可能收敛也可能不收敛.
对非线性方程组解的存在性的研究远不如线性方程组那样成熟，现有的解法也不象线性方程组那样有效.除极特殊的方程外，一般不能用直接方法求得精确解，目前主要采用迭代法求近似解.其中最著名的是牛顿法—拉卜森迭代法，也就是常说的牛顿法.而二分法的代码为：
erfen.m
function y=erfen(fun,a,b,esp)
if nargin<4 eps=1e-4;end
if feval_r(fun,a)*feval_r(fun,b)<0
n=1;
c=(a+b)/2;
while c>esp
if feval_r(fun,a)*feval_r(fun,c)<0
b=c;c=(a+b)/2;
elseif feval_r(fun,c)*feval_r(fun,b)<0
a=c;c=(a+b)/2;
else y=c;esp=10000;
end
n=n+1;
end
y=c;
elseif feval_r(fun,a)==0
y=a;
elseif feval_r(fun,b)==0
y=b;
else disp(‘these,may not be a root in the intercal’);
end
n
例3用二分法求方程x^2-x-1=0的正根，要求误差小于0.05.
解：首先编制函数文件：
Fc.m
function y=fc(x)
y=x^2-x-1;
return
在命令窗中输入：
>>erfen(‘fc’,0,10,0.05)
n =
56
ans =
1.6180
二分法在区间很大时，收敛接近根的速度很快，而当区间较小时，靠近速度变得相当缓慢，且计算量较大，对高精度要求的满足比较困难，因此，可以为其他各种迭代法提供迭代初值.
3 、牛顿法
牛顿法基本思想是将非线性问题逐步线性化而形成如下迭代程序：
（4）其中
是 (1)式的雅可比矩阵，x(0) 是方程(1)的解的初始近似.
这个程序至少具有2阶收敛速度.由x(k)算到x(k+1)的步骤为:
① 由x(k)算出f(x(k))及 ;
② 用直接法求线性方程组 的解△x(k)；
③ 求 .
由此看到迭代一次需计算1个分量函数值和1个分量偏导数值，并求解一次k阶线性方程组.为了避免切线法计算导数的麻烦，为此我们引进割线法.
由公式（4）得，为此可以用matlab来解非线性方程组的数值解法，其代码为：
function y=newtonpro(x0)
f1=f1ops;
x1=x0-fc(x0)/df(x0);
n=1;
while (norm(x1-x0)>1.0e-6)&(n<=100000000)
x0=x1;
x1=x0-fc(x0)/df(x0);n=n+1;
end
f2=flops;
f=f2-f1
x1
n
其中fc为原函数，df为导函数.
例1：构造不动点迭代求解如下方程组.
x1-0.7sin(x1)-0.2cos(x2)=0
x2-0.7cos(x1)+0.2sin(x2)=0
在(0.5,0.5)附近的解，迭代至1e-3的精度.
解：编制函数文件fc.m和导数文件df.m：
fc.m
function y=fc(x)
y(1)=x(1)-0.7*sin(x(1))-0.2*cos(x(2));
y(2)=x(2)-0.7*cos(x(1))+0.2*sin(x(2));
y=[y(1) y(2)];
return
df.m
function y=df(x)
y=[1-0.7*cos(x(1)) 0.2*sin(x(2))
0.7*sin(x(1)) 1+0.2*cos(x(2))];
return
于是在MATLAB命令窗中计算得：
>>x0=[0.5 0.5];
>>newtonpro(x0)
f =
875
x1 =
0.5265 0.5079
n =
12
可见Newton法收敛得要快一些，且可得到不同的解.但当函数有重根的过程相当缓慢.因此，当方程有重根时，要改用重根条件下的改进型Newton迭代法.
4、
割线法
上述牛顿法虽然收敛速度很快，但要求偏导数，这往往难以做到.而割线法虽采用了类似的思路但避免了求偏导数.这种方法的原理用一维方程最易说明.假定有一个一维问题f(x)=0，为了解此方程，有两个解的初始近似值x(k和x(k+1)，与之对应的函数值为f(k)=f[x(k)]和f(k-1)=f[x(k-1)].为了得到进一步的近似解x(k+1)，可以在点(x(k),f(k))和点(x(k-1),f(k-1))之间作一条直线，而这一直线与轴的交点就是下一个x(k+1).这可以考虑用差商(f(k)-f(k-1))/(x(k)-x(k-1))代替f(x(k))’，
迭代公式变为
x(k+1)=x(k)—f(k)* (x(k)-x(k-1))/ (f(k)-f(k-1))
（5）
将x(k+1)代入函数f(x)，可以产生出新的f(k+1)=f[x(k+1)]，于是，可以将最初的值x(k-1)舍弃，采用两个最新的x近似值继续迭代.
将以上构思推广到多维情况解非线性方程组，为减少解线性方程组次数，于是得到 称为解非线性方程组的割线法.其代码为：
ger.m
function y=ger(x0,x1)
x2=x1-fc(x1)*(x1-x0)/(fc(x1)-fc(x0));
n=1;
while (abc(x1-x0)>=1.0e-4)&(n<=100000000)
x0=x1;x1=x2;
x2=x1-fc(x1)*(x1-x0)/(fc(x1)-fc(x0));
n=n+1;
end
x2
n
return
例2：用割线法求方程f(x)=x^3-3*x-1在x0=2附近的根.误差限为1e-4，取x0=2，x1=1.9.
解：编制函数文件：
fc.m
function y=fc(x)
y=x^3-3*x-1;
>>ger(2,1.9)
x2=
1.8794
n =
4
从计算中可看出，割线法的收敛也是相当快的.理论证明可知割线法具有超线性的收敛性.
5、拟牛顿法
为减少牛顿法的计算量，避免计算雅可比矩阵及其逆，60年代中期出现了一类称为拟牛顿法的新算法,它有不同的形式，常用的一类是秩1的拟牛顿法，其中不求逆的程序为
式中
， ， ，称为逆拟牛顿公式.计算时先给出x(k)及 (k)
，由(9)逐步迭代到满足精度要求为止.每步只算1个分量函数值及2的计算量,比牛顿法一步计算量少得多.理论上已证明，当 x(k)及
(k) 选得合适时,它具有超线性收敛速度,但实践表明效率并不高于牛顿法.
其代码如下：
broyden.m
function y=broyden(x0)
%y=broyden(x0) x0为初值.
a=eye(length(x0));
x1=x0-fc(x0)/a;
n=1;
while
(norm(x1-x0)>=1.0e-6)&(n<=100000000)
x0=x1;
x1=x0-fc(x0)/a;
p=x1-x0;q=fc(x1)-fc(x0);
a=a+(q-p*a)’*p/norm(p);
n=n+1;
end
y=x1;
n
再解例1，其函数fc.m函数不变.
>>broyden(x0)
n =
22
ans=
0.5265 0.5079
由此例可用看出虽然拟牛顿法比Newton法收敛得慢一些，但在此方法中不必计算原函数的代数值，甚至连导数都不用，很方便.
6、不动点迭代法
设一元函数f是连续的，要解得方程是（1），为了进行迭代，变换方程形式为x=G(x)，于是构造迭代公式x(k+1)=G(x(k)).若当k趋于无穷大时,则x(k)趋于x*，则称此迭代为不动点迭代.易知构造不动点迭代的迭代公式不一定收敛速度也受所构造的迭代公式的好坏的影响.其函数如下：
iteratepro.m
function y=iteratepro(x)
x1=g(x);
n=1;
while (norm(x1-x)>=1.0e-6)&(n<=1000)
x=x1;
x1=g(x);n=n+1;
end
x1
n
return
对例1进行迭代得编制函数文件：
g.m
function y=g(x)
y(1)=x(1)-0.7*sin(x(1))-0.2*cos(x(2));
y(2)=x(2)-0.7*cos(x(1))+0.2*sin(x(2));
return
计算得：
>>iteratepro([0.5 0.5])
x1 =
0.5265 0.5079
n=
23
二 、总结
20世纪60年代中期以后，发展了两种求解非线性方程组(1)的新方法
一种称为区间迭代法或称区间牛顿法，它用区间变量代替点变量进行区间迭代，每迭代一步都可判断在所给区间解的存在惟一性或者是无解.这是区间迭代法的主要优点，其缺点是计算量大.另一种方法称为不动点算法或称单纯形法，它对求解域进行单纯形剖分，对剖分的顶点给一种恰当标号，并用一种有规则的搜索方法找到全标号单纯形,从而得到方程(1)的近似解.这种方法优点是，不要求f(
)的导数存在,也不用求逆，且具有大范围收敛性，缺点是计算量大.通过这几种解法，我们或多或少对求数值解法有一定的了解.至于以后对于非线性方程组的解法，要看哪个解法最有效最有用.
