2.3 数值分析
数值分析是计算数学的一个重要部分,它研究用计算机求解各种数学问题的数值计算方法,以及其理论与软件实现。用计算机求解数学问题通常经历以下5个步骤:
实际问题→数学模型→数值计算方法→程序设计→上机计算求出结果。
根据实际问题建立数学模型往往是应用数学的任务,计算数学更关注的是如何给出数值计算方法,并根据计算方法编制算法程序,从而求得最终的计算结果。
计算机及科学技术的快速发展使求解各种数学问题的数值方法也越来越多,解决问题的速度和效率也得到了很大的提升。数值分析涉计数学的各个分支,所包含的内容十分广泛,本节只介绍基本的数值计算方法及其理论内容,包括:插值与数据逼近;数值微分与积分;线性方程组的数值求解;非线性方程与方程组求解等。它以数学问题为导向,将理论与计算相结合,重点关注数学问题的数值算法及其理论,是一门依赖于计算机技术的实用性课程。
2.3.1 插值法与数值逼近
插值法有着悠久的历史,它来自古人的生产实践活动,可以追溯到一千多年前的隋唐时期,人们利用二次插值法制定了历法,进行天文计算。17世纪,微积分的产生给插值理论的建立提供了条件。近半个世纪以来,由于计算机的推广,插值法逐步延伸到造船、航空航天、机械加工等实际工程问题中,获得了更广泛的应用。常见的插值法有拉格朗日插值法、牛顿的等距节点插值法,以及均差插值公式、埃尔米特插值法、分段低次插值法、三次样条插值法等。
(1)拉格朗日插值。
若n次多项式lj(x)(j=0,1,…,n)在n+1个节点x0<x1<…<xn上满足式(2.1):
我们称这n+1个n次多项式l0(x),l1(x),…,ln(x)为节点x0,x1,…,xn上的n次插值基函数,可以表示为式(2.2):
则n次插值多项式可以表示为式(2.3):
由lk(x)的定义可得式(2.4):
此时的插值多项式Ln(x)称为拉格朗日插值多项式。
(2)均差与牛顿插值多项式。
利用插值基函数可以得到拉格朗日插值多项式,但当插值节点数发生增减变化时,需要重新进行计算。为了计算方便,可重新设计一种逐次生成插值多项式的方法。
已知f在插值点xi(i=0,1,…,n)上的值为f(xi)(i=0,1,…,n),要求n次插值多项式Pn(x)满足条件式(2.5):
则Pn(x)可以表示为式(2.6):
式中,a0,a1,…,an为待定系数。与拉格朗日插值不同的是,这里的Pn(x)是由基函数{1,x-x0,…,(x-x0)…(x-xn-1)}逐次递推得到的。为了确定a0,a1,…,an的表达式,需要引入均差的定义。
一般地,称式
为f(x)的k阶均差。
借助均差的定义,可以将x看作[a,b]上的一点,可得式(2.8)。
式中,
其系数ak满足式(2.9):
我们称Pn(x)为牛顿均差插值多项式,它比拉格朗日插值计算量小,更便于程序设计。
在数值计算中经常要计算函数值,当函数只在有限点集上给定函数值,要在包含该点集的区间上用公式给出函数的简单表达式,这就涉及在区间[a,b]上用简单函数逼近已知复杂函数的问题,这就是函数逼近问题。
2.3.2 数值积分与数值微分
积分和微分是用极限来定义的两种分析运算的,处理数值分析和数值微分的基本方法是逼近法:设构造某个简单函数P(x)近似f(x),然后对P(x)求积(求导)得到f(x)的积分(导数)近似值。插值求积公式分为牛顿-柯特斯公式和高斯公式两类,实际计算宜采用复合求积的方法。高斯公式精度高、计算稳定,但节点选取较困难。带权高斯求积方法能把复杂积分简单化,还可以直接计算奇异积分。基于理查森外推的龙贝格求积方法由于计算程序简单、精度较高,是一种在计算机上求积的有效算法。在数值微分中也有相似的算法,但由于计算的不稳定性,在数值微分中,步长的选取很重要。
2.3.3 解线性方程的直接方法与迭代法
关于线性方程组的数值解法一般有直接法和迭代法两类。直接法就是经过有限步算术运算,求得方程组的精确解,但由于实际计算中存在误差,这种方法有一定的局限性,只能求得线性方程组的近似解。迭代法是利用某种极限过程去逐步逼近线性方程组的精确解。迭代法要求计算机的存储单元较少、程序设计简单、原始系数矩阵在计算过程中始终不变,但存在收敛性及收敛速度的问题。迭代法是解决大型稀疏矩阵方程组的重要方法。
1.直接方法
(1)高斯消去法。
设有线性方程组:
或写为矩阵形式:
简记为Ax=b。
高斯消去法解方程组的基本思想,是用逐次消去未知数的方法把原方程组Ax=b化为与其等价的三角形方程组,而求解三角形方程组可用回代的方法求解。
下面我们讨论如何用算法来实现高斯消去的过程:
设A∈Rm×n(m>1),s=m:n(m-1,n),如果,利用高斯方法将A化为上梯形,且A(k)覆盖A,乘数mik覆盖aik。
对于k=1,2,…,s:
① 如果akk=0,则停止计算。
② 对于i=k+1,…,m
i.aik←mik=aik/akk
ii.对于j=k+1,…,n
aij←aij-mik*akj
由此可见,算法的第k步需要做m-k次除法和(m-k)(n-k)次乘法运算,因此,本算法大约需要s3/3-(m+n)s2/2+mns次乘法运算。当m=n时,总计大约需要n3/3次乘法运算。
(2)高斯主元素消去法
由高斯主元素消去法可知,在消元过程中可能出现的情况,即无法进行消去法,即使主元素,但很小时,用其作除数会导致其他元素数量级的严重增长和舍入误差的扩散,最后,计算结果不准确。
目前,主要使用的是列主元素消去法,并假定式(2.11)的A∈Rn×n为非奇异的。
设方程组(2.11)的增广矩阵为
首先,在A的第1列中选取绝对值最大的元素作为主元素,如;然后,交换B的第1行与第i1行,经第1次消元计算得。
重复上述过程,设完成第k-1步的选主元素,交换两行及消元计算,(A b)化为
式中,A(k)的元素仍记为aij,b(k)的元素仍记为bi。
第k步选主元素(在A(k)右下角方阵的第1列内选),即确定ik,
交换(A(k)b(k))第k行与ik(k=1,2,…,n-1)行的元素,再进行消元计算,最后将原线性方程组化为:
回代求解
下面,我们讨论如何用算法来实现列主元素消去法的过程。
设Ax=b,本算法用A具有行交换的列主元素消去法。消元结果冲掉A,乘数mij冲掉aij,计算解x冲掉常数项b,行列式存放在det中。
① det←1。
② 对于k=1,2,…,n-1,
i.按列选取主元素。
ii.如果aik,k=0,则停止计算(det(A)=0)。
iii.如果ik=k,则转iv,
akj↔aik,j(j=k,k+1,…,n)。
换行:bk↔bik,
det←det。
iv.消元计算。
对于i=k+1,…,n,
a.aik←mik=aik/akk。
b.对于j=k+1,…,n,
aij←aij-mik·akj。
c.bi←bi-mik·bk。
v.det←akk·det。
③ 如果ann=0,则停止计算[det(A)=0]。
④ 回代求解。
i.bn←bn/ann。
ii.对于i=n-1,…,2,1,
⑤ det←ann·det。
(2)迭代法。对于由工程技术中产生的大型稀疏矩阵方程组(A的阶数n很大,但零元素较多),利用迭代法求解线性方程组Ax=b更合适。下面主要介绍迭代法中的雅克比迭代法和高斯-塞德尔迭代法。
1)雅克比迭代法。
解Ax=b的雅克比迭代法的计算公式为
2)高斯-塞德尔迭代法。
解Ax=b的雅克比迭代法的计算公式为
下面我们讨论如何用算法来实现高斯-塞德尔迭代法的过程。
设Ax=b,其中A∈Rn×n为非奇异矩阵,且aii≠0(i=1,2,…,n),本算法用高斯-塞德尔迭代法解Ax=b,数组x(n)开始存放x(0),后存放x(k),N0为最大迭代次数。
① xi←0.0(i=1,2,…,n)。
② 对于k=1,2,…,N0,i=1,2,…,n,
迭代一次,这个算法需要的运算次数最多与矩阵A的非零元素的个数一样。