第三部分 常用数值计算方法

第三部分常用数值计算方法中国矿业大学环境与测绘学院1/13/2023测绘软件设计与实现内容概要数值分析研究对象与特点线性方程组解法曲线拟合与函数逼近一、数值分析研究对象与特点1.1数值分析的定义及其主要内容数值分析也称为计算方法，是计算数学的一个主要组成部分计算数学是数学科学的一个分支，主要研究用计算机求解各种数学问题的数值计算方法及其理论与软件实现数值分析的主要内容数值分析的内容包括函数的数值逼近、数值微分与数值积分、非线性方程数值解、数值线性代数、常微和偏微数值解法等1.1数值分析的定义及其主要内容虽然数值分析也是以数学问题为研究对象，但它不像纯数学那样只研究数学本身的理论，而是把理论与计算紧密结合，着重研究数学问题的数值方法及其理论数值分析不是各种数值方法的简单罗列和堆积，是一门内容丰富，研究方法深刻，有自身理论体系的课程数值分析既有纯数学高度抽象性与严密科学性的特点，又有应用数学的广泛性与实际试验的高度技术性的特点，是一门与计算机使用密切结合的实用性很强的数学课程1.2数值分析的特点面向计算机，能根据计算机特点提供切实可行的有效算法。有可靠的理论分析，能任意逼近并达到精度要求，对近似算法要保证收敛性和数值稳定性，还要对误差进行分析要有好的计算复杂性，时间复杂性好是指节省时间；空间复杂性好是指节省存储量，这也是建立算法要研究的问题，它关系到算法能否在计算机上实现要有数值实验，即任何一个算法除了从理论上要满足上述三点外，还要通过数值试验证明是行之有效的1.3数值计算的若干原则使用数值稳定的算法在运算过程中舍入误差不增长的算法称为数值稳定的算法，否则称为数值不稳定的算法避免绝对值很小的数作为除数避免两相近的数相减通过改变计算公式的形式或算法，可以避免或减少有效数字的损失1.3数值计算的若干原则防止大数吃掉小数编制程序时，要合理安排计算顺序，防止重要的参数被“吃掉”简化计算步骤，减少运算次数提高算法的运行效率二、线性方程组解法2.1引题已知线性方程组 a11x1+a12x2+…+a1nxn=b1

a21x1+a22x2+…+a2nxn=b2

……. an1x1+an2x2+…+annxn=bn求解x1、x2、…xn的值？2.2线性方程组的解法直接法在没有舍入误差的情况下，经有限四则运算求解其准确值的方法迭代法类似于方程求根的迭代法先给出一个解的初始近似值，按一定的法则逐步将其精化2.2线性方程组的解法直接法高斯消去法（主元消去法、高斯-约当消去法）矩阵分解法迭代法雅可比迭代法高斯-赛德尔迭代法2.3消去法高斯顺序消去法列主元素法高斯—约当消去法求矩阵的逆2.3.1高斯顺序消去法消去法就是按特定的顺序进行的矩阵初等变换方法，当消元按自然顺序进行时，称为高斯顺序消去法例：

2x1+x2+x3=7

4x1+5x2-x3=11

x1-2x2+x3=02.3.1高斯顺序消去法第一步：消元此线性方程组的增广矩阵为对增广矩阵进行初等行变换，使系数矩阵化为上三角形2.3.1高斯顺序消去法第二步：得同解方程组 2x1+x2+x3=7

3x2-3x3=-3

-2x3=-6第三步：对上三角形方程组进行回代求解

x3=(-6)/(-2)=3

x2=(-3+3x3)/3=2

x3=(7-x2-x3)/2=1高斯顺序消去法的计算机算法增广矩阵经k-1步消元后，增广矩阵化为高斯顺序消去法的计算机算法第k步消元：设，以第k行为基础，将以后各行中的化为0，为此先计算然后以第i行减去第k行乘以li,k，即于是得高斯顺序消去法的计算机算法经n-1步消元后，增广矩阵化为自下往上逐步回代即可求得真解高斯顺序消去法的计算机算法由行列式的初等变换与矩阵初等变换的关系可知，顺序消元法的每一步保持系数行列式不变，因此，原方程组之系数行列式的值为，可在求解过程中逐步累乘求得。高斯顺序消去法流程图2.3.2列主元素法单精度解方程组 10-9x1+x2=1

x1+x2=2/*精确解:

*/8个8个用高斯消去法计算：8个小主元可能导致计算失败。8个列主元素法可以有效解决这一问题！！2.3.2列主元素法高斯顺序消去法的弱点当主元时，消元过程不能继续进行；当主元时，消元过程可以进行，但是计算

时，出现很小的数作除数的现象，使得舍入误差增大，严重情况下导致解的失真。列主元素法为了提高算法的数值稳定性，在消去过程中，应选取绝对值尽可能大的元素作为主元，基于此思想实现的算法成为列主元素法。2.3.2列主元素法在消元的第k步，首先在第k列下放所有元素，

，…中取绝对值最大者为主元，设

则选为主元，若，则交换第r行与第k行，而后再进行消元。实际计算中，当很小时，求解结果将严重失真，继续计算已无意义，因此可设一阈值ε，当

时则停止计算。2.3.2.1列主元素法的基本步骤第一步：选主元：在ai1（i=1,…,n)中选取绝对值最大的元素作为主元素，即确定i，使交换行：即当i≠1时，交换第i行与第1行元素；消元计算：li1=ai1/a11（i=2,3,…，n）

aij←aij-lija11(i=2,3,…,n;j=i,i+1,…,n）第k步：选主元：在aik（i=k,…,n)中选取绝对值最大的元素作为主元素，即确定i，使交换行：即当i≠k时，交换第i行与第k行元素；消元计算：lkk=aik/akk（i=k,k+1,…，n）

aij←aij-lija11

（i=k,k+1,…,n;j=k,k+1,…,n）2.3.2.2列主元素法消去过程2.3.2.3严格对角占优矩阵严格对角占优矩阵满足如下条件：对角线上每一元素的绝对值均大于同行其他元素绝对值之和对角线上每一元素的绝对值均大于同列其他元素绝对值之和若矩阵A对称且严格对角占优，则消元过程中第k步主元必为，因此不必选主元，用顺序消去法即可2.3.3高斯—约当消去法在消元过程中选取主元后，先将主元化为1，而后将主元所在列上、下方各元素均化为0，这样消元的结果使系数矩阵化成了单位阵，无需回代就得到了原方程组的解，这种方法称之为高斯—约当消去法。 [A

|I

][I

|B]高斯—约当消去法（矩阵的初等行变换）2.4矩阵分解法—LU分解法高斯消元法的矩阵形式记记2.4矩阵分解法—LU分解法对比以上可见，消元过程相当于下述矩阵乘法运算

M(2)M(1)[A|b]=[U|y]由分块矩阵乘法可得 M(2)M(1)A=UM(2)M(1)b=y

令

L=(M(2)M(1))-1=(M(1))-1(M(2))-1由直接计算可知于是

A=LU，Ly=b2.4矩阵分解法—LU分解法可见只要消元过程能进行到底，就有以下等价关系据此，可以得到LU分解之后矩阵L和U的基本表达式2.4矩阵分解法—LU分解法例：求下述矩阵的LU分解所以2.4矩阵分解法—LU分解法有了矩阵A的LU分解计算公式，解线性方程组Ax=b就转化为依次解下三角方程组Ly=b与上三角方程组Ux=y，其计算公式如下：2.4矩阵分解法—LU分解法定理一

若矩阵A非奇异，则A能分解为LU的充分必要条件是A的顺序主子行列式不为0，即定理二

若非奇异矩阵A有LU分解，则此分解是唯一的。2.4.1LU分解实例已知

，

用LU分解法分别求解下列线性方程组

（1）Ax=b；

（2）A2z=b2.4.1LU分解实例解：

对方程组A2z=b，若先计算A2再求解，则计算量较大，在一般情况下，A2还会产生新的误差；

若令x=Az，则问题转化为依次解Ax=b，Az=x两个线性方程组；

因此，利用LU分解法，等价于求解四个方程组Ly=b，Ux=y，Lw=x，Uz=w，即在求出第（1）题之解x的基础上很快求出第（2）题之解z。2.4.1LU分解实例A的LU分解式解Ly=b，得解Ux=y，得第（1）题之解解Lw=x，得解Uz=w，得第（2）题之解作业&数值实验一1、试编程求解如下方程组（列主元素法）（1）（2）（3）

对（1）（2）两题观察每步消元结果的系数矩阵有何特点，右下方矩阵是否对称，列主元在何处等2、对第1题中的系数矩阵A进行LU分解，并用此分解法解对应的线性方程组2.5迭代法迭代法求解方程的根迭代法求方程组的解2.5.1迭代法求解方程的根迭代法的基本思想迭代法是一种重要的逐次逼近法将方程f(x)=0化为等价方程x=φ(x)，然后在隔根区间内取一点x0，按下式计算xk+1=φ(xk)(k=0,1,2…) （1）计算结果生成序列 x0,x1,…,xk,… （2）如果这个数列有极限，显然x*就是方程x=φ(x)之根。于是，可以从此数列中求得满足精度要求的近似根，这种求根方法称为迭代法，式（1）称为迭代格式，φ(x)称为迭代函数，x0称为迭代初值，数列（2）称为迭代序列如果迭代序列收敛，则称迭代格式（1）收敛，否则称迭代格式发散（1）

实例：迭代法求解方程的根例：用迭代法求方程x4+2x2-x-3=0在区间[1,1.2]内的实根解：对方程进行如下三种变形 x=φ1(x)=(3+x-2x2)1/4

x=φ2(x)=((x+4)1/2-1)1/2

x=φ3(x)=x4+2x2-3分别按以上三种形式建立迭代格式，并取x0=1进行迭代计算，结果如下xk+1

=φ1(xk)，x26=x27=1.124123xk+1

=φ2(xk)，x6=x7=1.124123xk+1

=φ3(xk)，x3=96,x4=8.495307×107准确根：x*=1.124123029，可见迭代格式不同，收敛情况也不同，第二种格式比第一种格式收敛快得多，第三种格式不收敛（2）收敛条件与误差估计定理一（收敛条件）设φ′(x)(在[a,b]上存在，且满足条件：

（1）当x∈[a,b]时，φ(x)∈[a,b]

（2）存在正数L<1，使对任意的x∈[a,b]，|φ′(x)|≤L<1

则

（1）方程x=φ(x)在[a,b]上有唯一根x*;

（2）对任意的迭代初值x0∈[a,b]，迭代序列xk+1=φ(xk)(k=0,1,2,…)收敛于x*。（2）收敛条件与误差估计定理一（收敛条件）设φʹ(x)(在[a,b]上存在，且满足条件：

（1）当x∈[a,b]时，φ(x)∈[a,b]

（2）存在正数L<1，使对任意的x∈[a,b]，|φʹ(x)|≤L<1

则

（1）方程x=φ(x)在[a,b]上有唯一根x*;

（2）对任意的迭代初值x0∈[a,b]，迭代序列xk+1=φ(xk)(k=0,1,2,…)收敛于x*。（2）收敛条件与误差估计定理一（收敛条件）设φʹ(x)(在[a,b]上存在，且满足条件：

（1）当x∈[a,b]时，φ(x)∈[a,b]

（2）存在正数L<1，使对任意的x∈[a,b]，|φʹ(x)|≤L<1

则

（1）方程x=φ(x)在[a,b]上有唯一根x*;

（2）对任意的迭代初值x0∈[a,b]，迭代序列xk+1=φ(xk)(k=0,1,2,…)收敛于x*。（2）收敛条件与误差估计定理二（误差预计）若在方程x=φ(x)之根x*的某邻域R={x||x-x*|<δ}内φʹ(x)存在，且存在正常数L<1，使

|φʹ(x)|≤L<1，

则任取x0∈R，迭代格式xk+1=(xk)均收敛于x*，且有下列误差估计式

反之，若在根x*的邻域R内|φʹ(x)|≥1

则迭代必发散（2）收敛条件与误差估计定理二给出了迭代格式在根x*邻近是否收敛的判定条件定理二强调，当条件|φʹ(x)|≤L<1成立时，迭代初值x0应取在根x*的邻域R中如果对于任意给定的x0，迭代格式均收敛，则称此格式具有全局收敛性对根x*的某邻域内的一点x0，迭代格式收敛，则称此格式具有局部收敛性在具体解题时，虽然无法判

别隔根区间是否是以x*为中

心的邻域，但只要他足够小，

在其中|φʹ(x)|≤L<1成立，则对

其中任取一点x0必能保证收敛（2）收敛条件与误差估计另外，定理2中给出的两个误差预计式也很有意义，由误差预计式可知，条件式中|φʹ(x)|≤L<1，L越小，迭代过程收敛越快；若L<1，但L≈1，则收敛很慢，此种格式不宜采用如：

例3中采用的三种迭代格式，在隔根区间(1,1.2)内，有

根据定理2可知，第三种格式发散，第1、2种迭代格式收敛，而第二种迭代格式比第一种迭代格式收敛要快很多，这与实际迭代结果完全吻合（2）收敛条件与误差估计根据误差预计式可知，若|xk-xk-1|<ε，则近似解xk有如下的误差预计因此，在正常收敛情况下，即L不十分接近于1时，可用|xk-xk-1|<ε控制迭代结束，这种用前后两次迭代结果估计误差的办法称为事后误差估计，实际应用中，控制迭代结束的条件也常取为

E<

ε

其中，2.5.2迭代法求方程组的解迭代法，由于它能充分避免系数矩阵中零元的存储和计算，因此，特别适用于解系数矩阵很高而非零元极少（即大型稀疏）的线性方程组首先将方程组Ax=b转化为等价方程组x=Mx+g，取初始向量x(0)，按下述格式迭代

x(k+1)=Mx(k)+g

生成向量序列x(1),x(2),

…,x(k),…，若此序列收敛于x*，则有x*=Mx*+g，即x*为原方程组之解。因此，可根据精度要求选择一个合适的x(k)作为近似解，其中，M称为迭代矩阵，若序列{x(k)}极限存在，则称此过程收敛，否则发散。（1）Jacobi迭代法设有n阶线性方程组 a11x1+a12x2+…+a1nxn=b1

a21x1+a22x2+…+a2nxn=b2 ……. an1x1+an2x2+…+annxn=bn其中，aii≠0(i=1,…,n)，则可将它改写为等价方程组（1）Jacobi迭代法建立迭代格式

或缩写为按此格式迭代求解的方法称为Jacobi迭代法，又称简单迭代法（1-1）Jacobi迭代法实例用Jacobi迭代法解线性方程组 10x1﹣x2-2x3=7.2 ﹣x1+10x2-2x3=8.3 ﹣x1﹣x2+5x3=4.2解：生成的Jaboci迭代格式

（1-1）Jacobi迭代法实例取x(0)=(0,0,0)T，迭代结果如下：

k

x1(k) x2(k) x3(k) 1 0.72 0.830001 0.84 2 0.971 1.07 1.15 …… 11 1.099993 1.199993 1.299991 12 1.099998 1.199998 1.299997方程的准确解：x*=(1.1,1.2,1.3)T，迭代结果误差不超过10-5（1）Jacobi迭代法为了引入Jacobi迭代格式的矩阵形式，将n阶线性方程组的系数分裂为A=D–L-U其中于是，Jacobi格式可写为矩阵格式x(k+1)=D-1(L+U)x(k)+D-1b其迭代矩阵MJ=D-1(L+U)（2）Gauss-Seidel迭代法Jacobi迭代法的缺点在求xi(k+1)时，前面已求出的最新分量x1(k+1),…,xi-1(k+1)未被利用直观情况下，应用迭代所得新值应比旧值更好Gauss-Seidel迭代法

（2）Gauss-Seidel迭代法Gauss-Seidel迭代可以缩写为将其写为矩阵格式

x(k+1)=D-1[Lx(k+1)+Ux(k)]+D-1b若把x(k+1)解出来，则得等价的迭代格式 x(k+1)=(D–L)-1Ux(k)+(D–L)-1b其迭代矩阵MG=(D-L)-1U（2-1）Gauss-Seidel迭代法实例用Gauss-Seidel迭代法解线性方程组 10x1﹣x2-2x3=7.2 ﹣x1+10x2-2x3=8.3 ﹣x1﹣x2+5x3=4.2解：生成的Gauss-Seidel迭代格式

（2-1）Gauss-Seidel迭代法实例取x(0)=(0,0,0)T，迭代结果如下：迭代8次，得x(8)=(1.099998,1.199999,1.3)T结论：用Gauss-Seidel迭代法比Jacobi迭代法收敛速度快注意：部分情况下，Jacobi迭代法收敛，Gauss-Seidel迭代法发散方程的准确解：x*=(1.1,1.2,1.3)T（3）

迭代法的收敛条件设x*是原方程组Ax=b的解，则它满足等价方程组

x*=Mx*+g (1)迭代序列{x(k)}满足

x(k+1)=Mx(k)+g (2)令

ε(k)=x(k)-x*则{x(k)}是否收敛于x*等价于向量序列{ε(k)}是否收敛于0(2)-(1)得

ε(k+1)=Mε(k)

（k=0,1,…）因此，ε(k)=Mε(k-1)

=M2ε(k-2)

=…=Mkε(0)因为x(0)是任取的，ε(0)=x(0)-x*是任意向量，因此ε(k)→0的充要条件是M(k)→0。（3）

迭代法的收敛条件根据上述分析，迭代格式收敛与否取决于迭代矩阵M的特性，经过严格证明可得如下基本定理定理

对任何初始向量x(0)及常向量g，迭代过程 x(k+1)=Mx(k)+g

收敛的充分必要条件是

其中，λi(i=1,2,…,n)为矩阵M的特征值，

ρ(M)称为M的谱半径（3-1）迭代法的收敛条件实例设方程组Ax=b的系数矩阵A为

试判别Jacobi迭代与Gauss-Seidel迭代是否收敛（3-1）迭代法的收敛条件实例解：Jacobi迭代矩阵

令f(λ)=det(λI-MJ)=0

f(λ)=

λ3-(27/16)λ+(27/32)=0

因为f(-1)>0,f(-2)<0，所以，f(λ)=0在(-2,-1)之间有一根，因此，，所以，Jacobi迭代不收敛！（3-1）迭代法的收敛条件实例Gauss-Seidel迭代矩阵为

令f(λ)=det(λI-MG)=0

f(λ)=

λ(λ2-(81/64)λ+(27/64))=0

得λ1=0，λ2,3=0.633±0.146i，进而，

因此，Gauss-Seidel迭代法收敛！数值实验二Jacobi迭代法与Gauss-Seidel迭代法的收敛性与收敛速度研究用Jacobi迭代法与Gauss-Seidel迭代法解下列Ax=b方程组的收敛性，通过上机计算，验证分析是否正确，并观察右端对迭代收敛是否有影响，比较两法的收敛速度（1）（2）（3）三、最小二乘曲线、曲面拟合3.1引题在科学实验和统计研究中，往往要从大量的实验数据(xi,yi)(i=0,1,…,m)中寻找其函数关系y=f(x)的近似表达式y=P(x)由于实验观测数据不可避免地带有误差，甚至是较大的误差，因此，在寻找y=P(x)的过程之中，要求其偏差ri=P(xi)-yi(i=0,1,…,m)按照某种标准最小，以反映所给数据的总体趋势，消除局部波动的影响，这就是曲线拟合问题，这样的函数P(x)称为拟合函数。3.2最小二乘法基本原理基于最小二乘法的曲线拟合问题的具体提法对给定的数据(xi,yi)(i=0,1,…,m)，在取定的函数类Φ中，求P(xi)∈Φ，使偏差ri=P(xi)-yi(i=0,1,…,m)的平方和为最小，即满足上式的函数P(x)叫最小二乘拟合函数或最小二乘解3.2.1多项式拟合对于给定的一组数据(xi,yi)(i=0,1,…,m)，求作n次多项式(n≤m)

使其满足

即将拟合函数取为多项式，这样的曲线拟合问题叫多项式拟合问题，Pn(x)称为最小二乘拟合多项式。3.2.1多项式拟合Q可以看作是a0,a1,…,an的多元函数，因此，上述拟合多项式的构造问题可以归结为多元函数的极值问题由多元函数极值的必要条件知，aj(j=0,1,…,n)满足即这是关于系数a0,a1,…,an的线性方程组3.2.1多项式拟合将上述关于系数的a0,a1,…,an线性方程组写成矩阵形式为称之为法方程组或正则方程组法方程组的一个明显特点是其系数矩阵为对称的，如果它非奇异，那么，法方程组必存在唯一的一组解。从中求出ak(k=0,1,…,n)，即可求得多项式拟合结果3.2.2一般最小二乘拟合实际应用中，针对所讨论问题的特点，拟合函数可能为其他类型如指数函数、三角函数、有理函数等，待定参数还可能出现在指数上、分母中等，这就属于一般最小二乘拟合问题（1）

线性最小二乘拟合设φ0(x),φ1(x),…,φn(x)为n+1个线性无关的函数，Φ为其所有线性组合生成的函数集合，记作Φ

=Span{φ0(x),φ1(x),…,φn(x)}，任取P(x)∈Φ，则有

P(x)关于参数a0,a1,…,an是线性的对给定的一组数据(xi,yi)(i=0,1,2,…,m)，在Φ中求P(x)，使其满足

这就是一般的线性最小二乘拟合问题（1）

线性最小二乘拟合同多项式拟合类似，上述问题归为多元函数极值问题，由多元函数极值的必要条件

即

它是关于参数a0,a1,…,an的线性方程组，写成