一种基于拟Newton法解非线性方程组的数值算法

1、一种基于拟newton法解非线性方程组的数值算法（吴志华 中国计量学院信息工程学院 310018 ）关键字：非线性方程组 牛顿法1、理论背景我们先考虑线性方程，线性方程组的解便不难得出了。与线性方程相比，非线性方程问题无论是从理论上还是从计算公式 上，都要复杂得多。对于一般的非线性方程= 计算方程的 根既无一定章程可寻也无直接法可言。例如，求解高次方程组 7f _兀3 +兀_1 5 = 0的根，求解含有指数和正弦函数的超越方程 k-cos（空）=0的零点。解非线性方程或方程组也是计算方法中的一 个主题。在解方程方面，牛顿（inewton）提出了方程求才艮的一种 迭代方法，被后人称为牛顿算法。三

2、百年来，人们一直用牛顿算法, 改善牛顿算法，不断推广算法的应用范围。牛顿算法，可以说是数值 计算方面的最有影响的计算方法。对于方程式/（兀）=°,如果/（兀）是线性函数，则它的求根是 容易的。牛顿法实质上是一种线性化方法，其基本思想是将非线 性方程式/（x）逐步归结为某种线性方程来求解。解非线性方程 组只是非线性方程的一种延伸和扩展。2.主要理论考虑方程组了(几£)= 0, 九("£) = °其中/,，九均为3,兀j多元函数。若用向量记号记兀=3,兀j w/r,f =(几,£)'，（°就可写成f(x) = 0.(2)

3、当，且加=1，,）中至少有一个是自变量兀（z=l，,/l）的非线性函数时，则称方程组为非线性方程 组。非线性方程组求根问题是前面介绍的方程即=1）求根的 直接推广，实际上只要把单变量函数/（兀）看成向量函数f（x） 则可将单变量方程求根方法推广到方程组（2）。若已给出方程组（2）的一个近似根 +" =（#,#），将函数尸（对的分量 /（兀）（i = 1,兀）在兀伙）用多元函数泰勒展开，并取其线性部分, 则可表示为f（x）« f（x）+）（x 兀）.令上式右端为零，得到线性方程组f'（x ）（x_x）=_f（兀伙）,其中dxn茁（兀）dx茁（兀）dx2dx（4） 称为

4、f（q为雅可比（jacobi）矩阵。求解线性方程组卩） 记解为兀"，则得卅+1）=卅）_f（*）tf（卅）伙=0丄）.（5） 这就是解非线性方程组（2）的牛顿法。3、算法牛顿法主要思想是用卅=x-f'（x）"f（x） 伙= 0,1,.）.进行迭代。因此首先需要算出foo的雅可比 矩阵fq）,再求过fq）求出它的逆当它达到某 个精度（x-幻时即停止迭代。具体算法如下:1. 首先宏定义方程组尸（兀），确定步长-和精度（x-k）；2. 求fs）的雅可比矩阵f'（x）；可用dxjxj + x ，xn） £ （兀，x-，xfj求出雅可比矩阵;3. 求雅可比矩

5、阵f3的逆f'' ；q .0、 将f（兀）右乘一个单位矩阵1°1丿，通过单位矩阵变 换实现求f©）的逆，用指针来存贮。4. 雅可比矩阵与其逆f'（x）|的相乘；5. 用q）来迭代；6. 当精度x-k大于x_k时，重复执行25步，直到精度小 于或等于x-k停止迭代，x_kj就是最后的迭代结果。其中 x_kl = j（铲）-护铲）胡）24、代码#include <iostreamh>#include <stdlib. h> #include <math> h> #include vconio. h> #de

6、f ine f 0 (xl, x2) (xl+2*x2-3)#define fl (xl,x2) (2*xl*xl+x2*x2-5)#define x_ 0.000001 #define matrixnum 2 int y=0;void main()int i, j,n;double p, *x;double *b;double *matrixf; 矩阵 fdouble *matrixf_;/矩阵f的雅可比矩阵的逆b= (double *)malloc(matrixnum):matrixf=(double *)malloc(matrixnum);matrixf_=(double *)mallo

7、c(matrixnum*matrixnum); cout<"请输入初值:"for (i=0; i<matrixnum; i+)cin»* (x+i);dop=0. 0;for (i=0; i<matrixnum; i+)*(b+i)=0;*matrixf=fo(*x, *(x+1);*(matrixf+1)=f1(*x, *(x+1); matrixf_=matrixf2(x);for (j=0; j<matrixnum; j+)* (b+i) +=* (matrixf_+ i*matrixnum+j) * (* (matrixf+j);*

8、 (x+i)=*(x+i)-*(b+i)；cout«* (x+i) «/z ”;cout«endl;for (i=0; i<matrixnum; i+)p+二pow(*(b+i), 2);y+;while(sqrt(p)>x_);cout«,z停止迭代，最终迭代结果为«*x«',' «* (x+l) «,/z«endl;delete matrixf;delete matrixf_;getcho ;double *matrixf2(double *x)int i, j；doubl

9、e t;double *matrixfl;/矩阵f的雅可比矩阵double *matrixf2;/矩阵f的雅可比矩阵的逆matrixfl= (double *) malloc (matrixnum*matrixnuin);matrixf2= (double *)malloc(matrixnum*matrixnum);for(i=0;i<matrixnum;i+)for (j=0;j<matrixnum;j+)if (i二二j)* (matr i xf2+i *mat r i xnum+j)=1;else *(matrixf2+i*matrixnum+j)=0;*matrixfl=(

10、fo(*x+x_), *(x+1)-f0(*x, *(x+1) )/x_;*(matrixfl+l) = (fo(*x, (*(x+1)+x_)-fo(*x, *(x+1)/ x_；*(matrixfl+2) = (f1(*x+x_), *(x+1)-fl(*x, *(x+l)/x_;* (matrixfl+3) = (f 1 (*x, (*(x+l)+x_) -fl (*x, *(x+l)/x_;/for(i=0;i<matrixnum;i+)/ cout*(x+i)<<endl;cout矩阵 f 在*x<,' «*(x+l)<的雅可比矩阵en

11、dl;for (i=0; i<matrixnum; i+)for (j=0; j<matrixnum; j+)cout* (matrixfl+i*matrixnum+j) ;cout<<endl;/求矩阵f的雅可比矩阵的逆t二*matrixfl;for (i=0, j=0;j<matrixnum;j+)*(matrixfl+i*matrixnum+j)/二t;*(matrixf2+i*matrixnum+j)/=t;t=*(matrixfl+l*matrixnum);for (i=l, j=0:j<matrixnum;j+)*(matr i xf1+i *m

12、atr i xnum+j)-二*(matrixfl+j)*t;*(matrixf2+i *matri xnum+j)-二*(matrixf2+j)*t;t=*(matrixfl+l*matrixnum+l);for (i=l, j=0;j<matrixnum;j+)* (matrixfl+i*matrixnum+j)/=t;*(matrixf2+i*matrixnum+j) /=t;for (i=i, j=0: j<matrixnum; j+)*(matrixfl+j)-二*(matrixfl+i*matrixnum+j)*t;*(matrixf2+j)-二*(matrixf2+i

13、*matrixnum+j)*t;for (i=0;i<matrixnum;i+)for (j=0;j<matrixnum;j+)cout*(matrixf1+i*matrixnum+j)" ； cout«endl;for (i=0;i<matrixnum;i+)for (j=0;j<matrixnum;j+)cout* (matrixf2+i*matrixnum+j)< " cout<<endl;c0ut«z/第 y« 次迭代结果为,' «* (x+1) «z/,«

14、endl;getcho ;return matrixf2;delete matrixf2;5、算法分析及改进措施牛顿法解非线性方程组是一种线性方法，即将非线性方程组 以一线性方程组来近似，由此构造一种迭代格式，用以逐次逼近 所求的解案。可以证明newton迭代至少是二阶收敛的，而且收敛速度快。 因此牛顿法是解非线性方程的常规方法。正因为newton法思想直观自然，是最常用的，也是研究其 它方法的出发点，该方法的不足恰好是其它方法研究的出发点。首先，newton法的每步迭代都要计算尸,它是由n 个偏导数值构造的矩阵，有些问题中每个值可能都很复杂，甚至 根本无法解析地计算。当n比较大时这部分是算法

15、中耗费时间 最多的，不仅如此，每步迭代还要解线性方程组fu）2-尸（无）当n很大时（如由离散非线性偏微方程导出的非线性方程组，可能有1°"13甚至更多），其工作量很大。其次，在许多情况下，初值兀0要有较严格的限制，在实际 应用中给出确保收敛的初值是十分困难的。非线性问题通常又是 多解的，给出收敛到所需要解的初值更加困难。再有,迭代过程中如果某一步xk处有fg 奇异或几乎奇异（后者指f（境）的条件数彳艮大），则newton法的计算将无法 进行下去。特别如果在尸（兀）=。的解兀处有尸（兀）奇异，不 仅计算困难，而且问题本身也变得十分复杂，以一无元代数方程 为例，这时方程产生重根

16、。为了克服newton法的上述缺点，我们可以采用newton法 和参数newton法克服前两种缺点，拟newton法可以克服第三种 缺点。这里就拟newton法为例叙述对牛顿法的改进k = 0,1,2, .我们用矩阵近似的代替广（耳），从而得到如下形式的迭代法: 母+1=忑-b力、（忑）其中均为非奇异的.为了不要每次迭代都计算逆矩阵，我们设法构造丹鸟直接逼近 广（兀）的逆矩阵广（兀）迭代公式化为：k = 0,1,2, .6. 根据实例分析结果对于如下非线性方程组£(兀1兀2)=坷+2兀2-3 = 0, f9 (%| %2)_ 5 = 0.用如上源程序运行。1.输入初值为1. 5 l

17、0进行迭代的结果分别为:第一次 1. 50. 75理论值分别为:1. 50. 75第二次 1.48810.755952第三次 1.48803 0.7559831.4880950. 7559521. 4880340. 755983停止迭代，最终迭代结果为1.4880340. 755983输入初值为2. 02. 0进行迭代的结果分别为:第一次第二次第三次第四次1.833331.527781.48871. 488030. 5833330.7361110. 7556520. 755983停止迭代，最终迭代结果为：1.488030.7559833.输入初值为1, 01. 5进行迭代的结果分别为：第一次1.90.55第二次1. 54220. 728901第三次1. 489250. 755376第四次1.488030.755983停止迭代，最终迭代结果为：1.488030.7559834.输入初值为1000010000进行迭代的结果分别为：第一次 9998.99-4997.99第二次 4999.66-2498. 33第十五次 1. 541740.729131第十六次 1.489230.753386第十七次 1. 488030. 755983停止迭代，最终迭代结果为：1.488030.755983以上