C语言解决多元多次方程
一 理论背景
我们先考虑线性方程，线性方程组的解便不难得出了。 与线性方程相比，非线性方程问题无论是从理论上还是从计算公式上，都要复杂得多。对于一般的非线性方程f (x ) =0，计算方程的根既无一定章程可寻也无直接法可言。例如，求解高次方程组7x -x +x -1.5=0的根，求解含有指数和正弦函数的超越
x
方程e -cos(πx ) =0的零点。解非线性方程或方程组也是计算方
63
法中的一个主题。在解方程方面，牛顿（I . Newton）提出了方程求根的一种迭代方法，被后人称为牛顿算法。三百年来，人们一直用牛顿算法，改善牛顿算法，不断推广算法的应用范围。牛顿算法，可以说是数值计算方面的最有影响的计算方法。
对于言程式
f (x ) =0, 如果f (x ) 是线性函数，则它的求根
是容易的。牛顿法实质上是一种线性化方法，其基本思想是将非线性方程式f (x ) 逐步归结为某种线性方程来求解。解非线性方程组只是非线性方程的一种延伸和扩展。 二 主要理论 考虑方程组
⎧f 1(x 1,... x n ) =0, ⎪
................. ⎨ (1) ⎪f (x ,... x ) =0.
n ⎩n 1
其中f 1,..., f n 均为(x 1,... x n ) 多元函数。若用向量记号记
x =(x 1,... x n ) T ∈R n , F =(f 1,..., f n ) T ，(1)
就可写成
F (x ) =0. (2)
当
n ≥2, , 且f i (i =1,..., n ) 中至少有一个是自变量
x i (i =1,..., n ) 的非线性函数时，则称方程组(1)为非线性方程
组。非线性方程组求根问题是前面介绍的方程即(n =1) 求根的直接推广，实际上只要把单变量函数f (x ) 看成向量函数F (x ) 则可将单变量方程求根方法推广到方程组(2)。若已给出方程组
k T
) , 将函数F (x ) 的分量(2)的一个近似根 x (k ) =(x 1k ,..., x n
f i (x )(i =1,..., n ) 在x (k ) 用多元函数泰勒展开，并取其线性部分，
则可表示为 F (x ) ≈F (x
(k )
) +F '(x (k ) )(x -x (k ) ).
令上式右端为零，得到线性方程组
(k ) (k ) (k )
'F (x )(x -x ) =-F (x ), (3)
其中
⎡∂f 1(x ) ∂f 1(x ) ⎢∂x ∂x 21 ⎢
⎢∂f 2(x ) ∂f 2(x ) ⎢∂x ∂x 21 ⎢F '(x ) =⎢⎢
⎢
⎢∂f n (x ) ∂f n (x ) ⎢
∂x 2⎣∂x 1
∂f 1(x ) ⎤
∂x n ⎥
⎥∂f 2(x ) ⎥∂x n ⎥⎥
⎥ (4) ⎥ ⎥∂f n (x ) ⎥
⎥∂x n ⎦
称为F (x ) 为雅可比（Jacobi ）矩阵。求解线性方程组(3)，并记解为x (k +1) ，则得 x
(k +1)
=x (k ) -F '(x (k ) ) -1F (x (k ) ) (k =0,1,...). (5)
这就是解非线性方程组(2)的牛顿法。
三．算法
(k +1)
=x (k ) -F '(x (k ) ) -1F (x (k ) ) 牛顿法主要思想是用x
(k =0,1,...). 进行迭代 。因此首先需要算出F (x ) 的雅可比矩阵
F '(x ) ，再求过F '(x ) 求出它的逆F '(x ) -1，当它达到某个精度
(x_k)时即停止迭代。
具体算法如下：
1． 首先宏定义方程组F (x ) ，确定步长x _和精度(x_k)； 2． 求F (x ) 的雅可比矩阵F '(x ) ；
可用
∂f i (x 1,..., x j ,..., x n )
∂x j
求出雅可比矩阵；
=
f i (x 1,..., x j +x _,...,x n ) -f i (x 1,..., x j ,..., x n )
x _
3． 求雅可比矩阵F '(x ) 的逆F '(x ) ；
⎛1 0⎫ ⎪
将F '(x ) 右乘一个单位矩阵 ⎪，通过单位矩阵变
0 1⎪⎝⎭
-1
换实现求F '(x ) 的逆，用指针来存贮。
4． 雅可比矩阵F '(x ) 与其逆F '(x ) 的相乘； 5． 用(5)来迭代；
6． 当精度x_ki 大于x_k时，重复执行2——5步，直到精度
小于或等于x_k停止迭代，x_ki 就是最后的迭代结果。
其中x_ki =
-1
四．代码
#include #include #include #include
#define f0(x1,x2) (x1+2*x2-3) #define f1(x1,x2) (2*x1*x1+x2*x2-5) #define x_ 0.000001 #define matrixNum 2
double *matrixF2(double *x); int y=0; void main() {
int i,j,n; double p,*x; double *b;
double *matrixF; //矩阵F
double *matrixF_; //矩阵F 的雅可比矩阵的逆 b=(double *)malloc(matrixNum);
matrixF=(double *)malloc(matrixNum);
matrixF_=(double *)malloc(matrixNum*matrixNum);
cout
for(i=0;i
cin>>*(x+i); do
{
p=0.0;
for(i=0;i
*(b+i)=0;
*matrixF=f0(*x,*(x+1));
*(matrixF+1)=f1(*x,*(x+1)); matrixF_=matrixF2(x);
for(i=0;i
for(j=0;j
*(b+i)+=*(matrixF_+i*matrixNum+j)*(*(matrixF+j)); *(x+i)=*(x+i)-*(b+i);
cout
cout
for(i=0;i
p+=pow(*(b+i),2);
y++;
}while(sqrt(p)>x_); cout
最
终
迭
代
结
果
为
"
double *matrixF2(double *x) {
int i,j; double t;
double *matrixF1; //矩阵F 的雅可比矩阵
double *matrixF2; //矩阵F 的雅可比矩阵的逆 matrixF1=(double *)malloc(matrixNum*matrixNum); matrixF2=(double *)malloc(matrixNum*matrixNum); for(i=0;i
if(i==j)
*(matrixF2+i*matrixNum+j)=1;
else *(matrixF2+i*matrixNum+j)=0;
*matrixF1=(f0((*x+x_),*(x+1))-f0(*x,*(x+1)))/x_;
*(matrixF1+1)=(f0(*x,(*(x+1)+x_))-f0(*x,*(x+1)))/x_;
*(matrixF1+2)=(f1((*x+x_),*(x+1))-f1(*x,*(x+1)))/x_;
*(matrixF1+3)=(f1(*x,(*(x+1)+x_))-f1(*x,*(x+1)))/x_;
//for(i=0;i
cout
for(i=0;i
for(j=0;j
cout
cout
}
//求矩阵F 的雅可比矩阵的逆
t=*matrixF1;
for(i=0,j=0;j
{
*(matrixF1+i*matrixNum+j)/=t;
*(matrixF2+i*matrixNum+j)/=t;
}
t=*(matrixF1+1*matrixNum);
for(i=1,j=0;j
*(matrixF1+i*matrixNum+j)-=*(matrixF1+j)*t;
*(matrixF2+i*matrixNum+j)-=*(matrixF2+j)*t;
}
t=*(matrixF1+1*matrixNum+1);
for(i=1,j=0;j
{
*(matrixF1+i*matrixNum+j)/=t;
*(matrixF2+i*matrixNum+j)/=t;
}
t=*(matrixF1+1);
for(i=i,j=0;j
*(matrixF1+j)-=*(matrixF1+i*matrixNum+j)*t;
*(matrixF2+j)-=*(matrixF2+i*matrixNum+j)*t;
}
for(i=0;i
for(j=0;j
cout
cout
}
for(i=0;i
for(j=0;j
cout
cout
cout
第
"
次
迭
代
结
果
"
getch();
return matrixF2; delete [] matrixF1; delete [] matrixF2; }
五．算法分析及改进措施
为
牛顿法解非线性方程组是一种线性方法，即将非线性方程组
以一线性方程组来近似，由此构造一种迭代格式，用以逐次逼近 所求的解案。
可以证明Newton 迭代至少是二阶收敛的，而且收敛速度快。
因此牛顿法是解非线性方程的常且方法。
正因为Newton 法思想直观自然，是最常用的，也是研究其
它方法的出发点，该方法的不足恰好是其它方法研究的出发点。 首先，Newton 法的每步迭代都要计算F '(x k ) ，它是由n
个偏导数值构造的矩阵，有些问题中每个值可能都很复杂，甚至 根本无法解析地计算。当n 比较大时这部分是算法中耗费时机 最多的，不仅如此，每步迭代还要解线性方程组
2
F '(x k ) x =-F (x k ) ，
当n 很大时（如由离散非线性偏微方程导出的非线性方程组，
n 可能有10
4
106甚至更多），其工作量很大。
其次，在许多情况下，初值x 0要有较严格的限制，在实际
应用中给出确保收敛的初值是十分困难的。非线性问题通常又是 多解的，给出收敛到所需要解的初值更加困难。
再有，迭代过程中如果某一步x k 处有F '(x k ) 奇异或几乎奇
异（后者指F '(x k ) 的条件数很大），则Newton 法的计算将无法
**
'F (x ) 奇异，不x F (x ) =0进行下去。特别如果在的解处有
仅计算困难，而且问题本身也变得十分复杂，以一无元代数方程为例，这时方程产生重根。
为了克服Newton 法的上述缺点，我们可以采用Newton 法
和参数Newton 法克服前两种缺点，拟Newton 法可以克服第三种缺点。
这里就拟Newton 法为例叙述对牛顿法的改进
我们用矩阵B k 近似的代替f '(x k ) ，从而得到如下形式的迭代法:
x k +1=x k -B k f (x k ) k =0,1,2, ….
其中B k 均为非奇异的.
为了不要每次迭代都计算逆矩阵，我们设法构造H k 直接逼近f '(x k ) 的逆矩阵f '(x k ) ，迭代公式化为 :
x k +1=x k -H k f (x k ) k =0,1,2, ….
六．根据实例分析结果
对于如下非线性方程组 -1-1
⎧⎪f 1(x 1, x 2) =x 1+2x 2-3=0,
⎨ 22f (x x ) =2x +x -5=0. ⎪12⎩21, 2
用如上源程序运行。
1．输入初值为1.5 1.0
进行迭代的结果分别为： 理论值分别为： 第一次 1.5 0.75 1.5 0.75 第二次 1.4881 0.755952 1.488095 0.755952 第三次 1.48803 0.755983 1.488034 0.755983 停止迭代，最终迭代结果为1.488034 0.755983
2. 输入初值为2.0 2.0
进行迭代的结果分别为：
第一次 1.83333 0.583333
第二次 1.52778 0.736111
第三次 1.4887 0.755652
第四次 1.48803 0.755983
以上表明
1）迭代中存在误差；
这是由于求雅可比式时用差商法来近似替代了，要过一系列的运算，从而误差在所难免。由于精度x_取0.000001，已经足够小了，所以迭代值与理论值相差并不大。
2）迭代次数与初值有关；
当初值与真实值越接近时迭代次数越少。
3）这种牛顿法迭代算法收敛性比较好，基本上能达到较强的收敛
要求。
七．参考文献
[1] 李庆扬，王能超，易大义 编 《数值分析》 清华大学
出版社，施普林格出版社 2001年8月
[2] 关治，陆金甫 编 《数值分析基础》 高等教育出版社 1998年5月
[3] 邓建中，葛仁杰，程正兴 编 《计算方法》 西安交通
大学出版社
[4] 王则柯 编 《计算的复杂性》 湖南教育出版社