第 16 讲数值计算

第16讲数值计算 • 掌握多项式的计算方法 • 掌握多元一次代数方程组的求根计算 • 掌握非线性方程组的求根计算 • 掌握积分计算的方法 • 掌握求逆矩阵的方法

计算机数据处理 数值数据处理（数值运算）：主要是数学计算，即求数值解，如求方程的根、求函数的定积分等非数值数据处理（又称非数值运算）：最常见的是事务管理领域，如图书馆管理、银行管理、自动化生产线管理等

多项式的计算 设n次多项式的通用形式为： Pn(x) = a0 + a1x + a2x2 + …… + anxn (其中n≥1) 给定一个x的值，求多项式函数Pn(x)的值。注意：在使用计算机运行一个算法时，乘除法要比加减法消耗时间长，因此在编写算法时，总是试图减少乘除法的次数。

1 逐项递推算法 • 逐项递推计算方法可以归结为如下关系式： • tk = x*tk-1 k=0，1，2，3 ……n • uk = uk-1 + ak*tk • t0 = 1 • u0 = a0 • 该算法避免了x的k次幂重复计算，从而提高了整个算法的效率

double Compute_Poly _1( double a[], int n , double x ) { double result , t ; int i ; t = x ; result = a[0] + a[1] * t ; for ( i = 2 ; i < = n ; i ++ ) { t = t * x ; result = result + a[i] * t ; } return result ; }

2 秦九韶算法 将多项式按降次幂写成如下形式： Pn(x) = (……((anx+an–1)x +……+a1)x+a0 仔细观察，不难发现一种计算方法，从里向外一层一层地计算，即有如下递推关系式： u k = a n（k = n –1，……，2，1，0） u k = u k + 1 * x + a k

double Compute_Poly _2( double a[], int n , double x ) { double result ; int i ; result = a[n] ; for ( i = n - 1 ; i > = 0 ; i–– ) { result = result * x + a[i] ; } return result ; } 秦九韶算法的运算效率较高，秦九韶算法也称为霍纳（Horner）方法。

多元一次代数方程组的求根计算 含多个未知元的线性方程组的求解问题，即求 a11 x1 + a12 x2 + …… + a1n xn = b1 a21 x1 + a22 x2 + …… + a2n xn = b2 …………………………………… an1 x1 + an2 x2 + …… + ann xn = bn 的解x1,x2,x3,…,xn的值，其中aij和bi均为常数。将系数、未知元和bi写成矩阵形式： AX = B 当detA≠0时，多元一次方程组的解存在且唯一。这里介绍约当消元法和迭代法求解多元一次方程组。

1 约当消元法 • 约当消元法的主要步骤： • 选主元 • 列交换 • 行交换 • 归一化 • 消元化

int Djordan(double a[][n],double b[]) { int i,j,k,l,flag=1, js[n]; //数组js[]用于记录列交换信息 double d,temp; for(i=0;i<=n-1;i++) js[i]=i; //先记录原方程的列信息 for(k=0;k<=n-1;k++) { d=0.0; for(i=k;i<=n-1;i++) for(j=k;j<=n-1;j++) if(absd(a[i][j]>d)) { d=absd(a[i][j]); js[k]=j; l=i; }; if(d == 0.0) { flag=0; return flag; //方程无解，又称奇异 } if(js[k]!=k) //列交换 for(i=0;i<=n-1;i++) { temp=a[i][k]; a[i][k]=a[i][js[k]]; a[i][js[k]]=temp; }

if( l!=k ) //行交换 { for(j=k;j<=n-1;j++) { temp=a[k][j]; a[k][j]=a[l][j]; a[l][j]=temp; } temp=b[k]; b[k]=b[l]; b[l]=temp; } for(j=k+1;j<=n-1;j++) a[k][j]/=a[k][k]; //归一化 b[k]/=a[k][k]; for(i=0;i<=n-1;i++) //消元化 { if(i!=k) { for(j=k+1;j<=n-1;j++) a[i][j]-=a[i][k]*a[k][j]; b[i]-=a[i][k]*b[k]; } } } for(k=n-2;k>=0;k--) if (js[k]!=k) { temp=b[k]; b[k]=b[js[k]]; b[js[k]]=temp; } flag=1; return flag; }

2 迭代法 迭代法是指：把多元一次方程组的解看作是某种极限过程的极限，而在实现这一极限过程每一步的结果是把前一步所得的结果施行相同的演算步骤得到。注意：用迭代方法时，不可能将极限过程进行到底，而只能把迭代进行有限多次

迭代法步骤 将方程组做如下简单变形： a11 x1 = b1 ― ( + a12 x2 + a13x3 + …… + a1n xn ) a22 x2 = b2 ― ( a21 x1 + a23x3 + …… + a2n xn) ………………………………………………………… ann xn = bn ― ( an1 x1 + an2 x2 + a23x3 + …… ) （假设aii都不为0 ）

塞德尔（Seidel）迭代算法步骤及算法： 第1步：任取一组初值作为方程的根；第2步：当计算精度不够，并且迭代次数大于0时，循环做第3和第4步，否则退出循环，转到第5步；第3步：将一组根Xm―1代入方程组计算得到Xm；第4步：迭代次数减1；第5步：判断是否计算出方程的根。

int Seidel_compute ( double a[][n], double x[] , double b[] , double e ) // a为系数矩阵，n为未知元的个数，e为计算精度，b[]为常数向量，x[]为根向量 { double p,s,t; int i,j,L=500; //假设迭代500次 for(i=0;i<=n-1;i++) x[i]=0; p=e+1; while ((p>=e)&&(L>0)) { p=0.0; for(i=0;i<=n-1;i++) { t=x[i]; s=0; for(j=0;j<=n-1;j++) //等号右边的计算 if(j!=i) { s+=a[i][j]*x[j]; };

x[i]=(b[i]-s)/ a[i][i]; //迭代解出新的xi if(absd(x[i]-t)>p) p=absd(x[i]-t); } L--; } if((p>=e)&&(L==0)) return 0; // 方程无解 else return 1; // 方程有解 }

求逆矩阵 通过矩阵运算来求解多元一次方程组：设有非奇异矩阵A，如果有下式成立，则称A―1为逆矩阵： AA-1 = A-1A = E 设A为方程等号左边系数矩阵，B为右边常数向量，X为所有未知变元的向量。则有 A X = B A―1 A X = A―1 B （A―1 A）X = A―1 B X = A―1 B

求逆矩阵的主要步骤如下： 第1步：在右下方阵中选主元；第2步：如果主元不是a kk做行列交换；第3步：a kk = 1 / a kk；第4步：a kj = a kj * a kk( j = 1，2，3，…，n；j ≠ k )；第5步：a ij = a ij ― a ik * a kj(i = 1，2，…，n,i≠k； j = 1，2，…，n，j≠k )；第6步：a ik = ― a ik * a kk( i = 1，2，3，…，n；i ≠ k ) 第7步：恢复原方程行列次序。选主元的意义在于：对于某个k，若a kk = 0，则运算无法进行下去。而当| a kk |很小时，计算过程受精度影响不稳定

#include<iostream.h> #define n 4 int dinv(double a[][n]) //全选主元矩阵求逆函数 { int js[n],is[n]; double temp; int i,j,k; double d; for(i=0;i<=n-1;i++) { js[i]=i; is[i]=i; } for(k=0;k<=n-1;k++) { d=0.0; for(i=k;i<=n-1;i++) for(j=k;j<=n-1;j++) if(absd(a[i][j])>d) { d=absd(a[i][j]); is[k]=i; js[k]=j; };

if(d==0.0) { cout<<"奇异\n"; return 1; //奇异 } if(is[k]!=k) //行交换 for(j=0;j<=n-1;j++) { temp=a[k][j]; a[k][j]=a[is[k]][j]; a[is[k]][j]=temp; } if(js[k]!=k) //列交换 for(i=0;i<=n-1;i++) { temp=a[i][k]; a[i][k]=a[i][js[k]]; a[i][js[k]]=temp; } a[k][k]=1/a[k][k]; for(j=0;j<=n-1;j++) if (j!=k) a[k][j]*=a[k][k]; for(i=0;i<=n-1;i++) if(i!=k) { for(j=0;j<=n-1;j++) if(j!=k) a[i][j]-=a[i][k]*a[k][j]; };

for(i=0;i<=n-1;i++) if(i!=k) { a[i][k]=-a[i][k]*a[k][k]; }; } for(k=n-1;k>=0;k--) { for(j=0;j<=n-1;j++) //行恢复 if (js[k]!=k) { temp=a[k][j]; a[k][j]=a[js[k]][j]; a[js[k]][j]=temp; } for(i=0;i<=n-1;i++) //列恢复 if (is[k]!=k) { temp=a[i][k]; a[i][k]=a[i][is[k]]; a[i][is[k]]=temp; } } return 0; }

测试主函数 • int main() • { double a[n][n]={1,2,1,-2,2,5,3,-2,-2,-2,3,5,1,3,2,3}; • double b[n]={4,7,-1,0}; • double x[n]; • int i,k; • // double a[n][n]={20,2,3,1,8,1,2,-3,15}; • // double b[n]={24,12,30}; • dinv(a); • for(i=0;i<=n-1;i++) • { • x[i]=0.0; • for(k=0;k<=n-1;k++) x[i]+=a[i][k]*b[k]; • } • for(i=0;i<=n-1;i++) cout<<"\t"<<x[i]; • cout<<endl; • return 0; • }

积分计算 求积分算法的基本思想是：将积分区间N等份，分别计算梯形面积，并规定计算精度，N的大小取决于计算精度。 N取值越大，计算的积分值越精确。

#include<iostream.h> #include<stdio.h> double f(double x); double hsab(double a,double b,double E) { int n,k; double h; double s,t; n=1; h=b-a; s=h*(f(a)+f(b))/2; t=s+E+1; while (absd(s-t)>=E) { t=s; s=0; for(k=0;k<=n-1;k++) s=s+f(a+k*h+0.5*h); s=t/2+h*s/2; h/=2; n*=2; } return s; }

double f(double x) { double result; result = 4 / (1 + x*x); return result; } int main() { cout<<"S = "<<hsab(0,1,0.0000001)<<endl; return 0; } 测试主函数

作业1：本章练习题1 作业2：本章练习题2 作业2：本章练习题4 作业3：本章练习题6

第 16 讲 数值计算