tceic.com
学霸学习网 这下你爽了
赞助商链接
当前位置:首页 >> 高等教育 >>

数值计算方法


数 值 计 算 方 法

第 1 章 绪 论.......................................................................................................................................................... 4 1.1 数值计算方法的对象与特点.................................................................................................................. 4 1.1.1 1.1.2 1.2 什么是数值计算方法.................................................................................................................. 4 数值计算方法的内容.................................................................................................................. 4

误差来源与误差分析.............................................................................................................................. 5 1.2.1 1.2.2 误差的来源.................................................................................................................................. 5 绝对误差、相对误差与有效数字...............................................................................................6

1.3

误差传播与若干防治办法.................................................................................................................... 11

第 2 章 非线性方程求根...................................................................................................................................... 14 2.1 实根的对分法........................................................................................................................................ 14 2.1.1 2.1.2 2.2 2.3 2.4 2.5 2.6 2.7 使用对分法的条件.................................................................................................................... 14 对分法求根算法........................................................................................................................ 15

迭代法.................................................................................................................................................... 16 迭代法的加工........................................................................................................................................ 20 牛顿迭代法............................................................................................................................................ 21 弦截法.................................................................................................................................................... 25 非线性方程组的牛顿方法.................................................................................................................... 27 程序示例................................................................................................................................................ 31

第 3 章 解线性方程组的直接法.......................................................................................................................... 35 3.1 消元法.................................................................................................................................................... 36 3.1.1 3.1.2 三角形方程组的解.................................................................................................................... 36 高斯消元法与列主元消元法.....................................................................................................38

3.1.3 高斯-若尔当(Gauss-Jordan)消元法................................................................................. 47 3.2 直接三角分解法.................................................................................................................................... 49 3.2.1 3.2.2 3.2.3 3.3 范 3.3.1 3.3.2 3.4 多利特尔分解............................................................................................................................ 52 追赶法........................................................................................................................................ 58
T 对称矩阵的 LDL 分解................................................................................................................ 59

数.................................................................................................................................................... 63 向量范数.................................................................................................................................... 63 矩阵范数.................................................................................................................................... 66

矩阵的条件数........................................................................................................................................ 71

3.5 3.6

病态方程组的解法................................................................................................................................ 72 程序示例................................................................................................................................................ 73

第 4 章 解线性方程组的迭代法.......................................................................................................................... 87 4.1 雅可比(Jacobi)迭代法.................................................................................................................... 88 4.1.1 4.1.2 4.2 4.3 4.4 4.5 雅可比迭代格式........................................................................................................................ 89 雅可比迭代收敛条件................................................................................................................ 92

高斯-塞德尔(Gauss-Seidel)迭代法...............................................................................................94 逐次超松弛(SOR)迭代法.................................................................................................................. 99 逆矩阵计算.......................................................................................................................................... 103 程序示例.............................................................................................................................................. 104 值.................................................................................................................................................... 113

第5章 插 5.1 5.2

引言...................................................................................................................................................... 113 拉格朗日(Lagrange)插值.............................................................................................................. 114 5.2.1 5.2.2 5.2.3 线性插值.................................................................................................................................. 114 二次插值.................................................................................................................................. 116 n 次拉格朗日插值多项式....................................................................................................... 120

5.3 5.4 5.5

埃特金逐步插值法.............................................................................................................................. 126 埃尔米特(Hermite)插值................................................................................................................ 130 分段插值.............................................................................................................................................. 138 5.5.1 5.5.2 龙格(Runge)现象................................................................................................................ 138 分段线性插值.......................................................................................................................... 139

5.6

三次样条函数...................................................................................................................................... 141 5.6.1 5.6.2 三次样条插值的 M 关系式...................................................................................................... 141 三次样条插值的 m 关系式...................................................................................................... 145

5.7

程序示例.............................................................................................................................................. 147

第 6 章 曲线拟合的最小二乘法........................................................................................................................ 155 6.1 6.2 6.3 6.4 6.5 拟合曲线.............................................................................................................................................. 155 线性拟合和二次拟合函数.................................................................................................................. 157 解矛盾方程组...................................................................................................................................... 161 权.......................................................................................................................................................... 169 用正交多项式作最小二乘拟合.......................................................................................................... 171

6.6

程序示例.............................................................................................................................................. 176

第 7 章 数值微分和数值积分............................................................................................................................ 183 7.1 数值微分.............................................................................................................................................. 183 7.1.1 7.1.2 7.1.3 7.2 差商与数值微分...................................................................................................................... 183 插值型数值微分...................................................................................................................... 187 样条插值数值微分.................................................................................................................. 189

数值积分.............................................................................................................................................. 189 7.2.1 7.2.2 7.2.3 插值型数值积分...................................................................................................................... 190 牛顿-柯特斯(Newton-Cote’s)积分.................................................................................192 求积公式的收敛性与稳定性...................................................................................................198

7.3

复化数值积分...................................................................................................................................... 199 7.3.1 7.3.2 7.3.3 7.3.4 复化梯形积分.......................................................................................................................... 199 复化辛普森积分...................................................................................................................... 201 复化积分的自动控制误差算法...............................................................................................203 龙贝格(Romberg)积分........................................................................................................ 207

7.4 7.5 7.6

重积分计算.......................................................................................................................................... 209 高斯(Gauss)型积分公式介绍........................................................................................................ 214 程序示例.............................................................................................................................................. 218

第 8 章 常微分方程数值解................................................................................................................................ 224 8.1 欧拉(Euler)公式............................................................................................................................ 225 8.1.1 8.1.2 8.1.3 8.2 基于差商的欧拉公式.............................................................................................................. 225 欧拉公式的收敛性.................................................................................................................. 229 基于数值积分的差分方法...................................................................................................... 232

龙格-库塔方法.................................................................................................................................... 234 8.2.1 二阶龙格-库塔方法................................................................................................................ 234

8.3 8.4

线性多步法.......................................................................................................................................... 241 常微分方程组的数值解法.................................................................................................................. 245 8.4.1 8.4.2 一阶常微分方程组的数值解法...............................................................................................245 高阶常微分方程数值方法...................................................................................................... 248

8.5 8.6

常微分方程的稳定性.......................................................................................................................... 249 程序示例.............................................................................................................................................. 252

第1章
1.1

绪 论

数值计算方法的对象与特点

1.1.1 什么是数值计算方法
现代的科学技术发展十分迅速,他们有一个共同的特点,就是都有大量的数据问题。 比如,发射一颗探测宇宙奥秘的卫星,从卫星设计开始到发射、回收为止,科学家和工程技术人员、 工人就要对卫星的总体、部件进行全面的设计和生产,要对选用的火箭进行设计和生产,这里面就有许许 多多的数据要进行准确的计算。发射和回收的时候,又有关于发射角度、轨道、遥控、回收下落角度等等 需要进行精确的计算。 有如,在高能加速器里进行高能物理试验,研究具有很高能量的基本粒子的性质、它们之间的相互作 用和转化规律,这里面也有大量的数据计算问题。 计算问题可以数是现代社会各个领域普遍存在的共同问题,工业、农业、交通运输、医疗卫生、文化 教育等等,各行各业都有许多数据需要计算,通过数据分析,以便掌握事物发展的规律。 研究计算问题的解决方法和有关数学理论问题的一门学科就叫做计算方法。计算方法属于应用数学的 范畴,它主要研究有关的数学和逻辑问题怎样由计算机加以有效解决。

1.1.2 数值计算方法的内容
数值计算方法也叫做计算数学或数值分析。数值计算方法主要内容包括非线性方程求根、线性代数方 程组解法、微分方程的数值解法、插值问题、函数的数值逼近问题、概率统计计算问题等等,还要研究解 的存在性、惟一性、收敛性和误差分析等理论问题。 我们知道五次及五次以上的代数方程不存在求根公式,因此,要求出五次以上的高次代数方程的解, 一般只能求它的近似解,求近似解的方法就是数值分析的方法。对于一般的超越方程,如对数方程、三角 方程等等也只能采用数值分析的办法。怎样找出比较简洁、误差比较小、花费时间比较少的计算方法是数 值分析的主要课题。 在求解方程的办法中,常用的办法之一是迭代法,也叫做逐次逼近法。迭代法的计算是比较简单的, 是比较容易进行的。 迭代法还可以用来求解线性方程组的解。 求方程组的近似解也要选择适当的迭代公式, 使得收敛速度快,近似误差小。 在线性代数方程组的解法中,常用的有塞德尔迭代法、共轭斜量法、超松弛迭代法等等。此外,一些 比较古老的普通消去法,如高斯法、追赶法等等,在利用计算机的条件下也可以得到广泛的应用。 在数值计算方法中,数值逼近也是常用的基本方法。数值逼近也叫近似代替,就是用简单的函数去代 替比较复杂的函数,或者代替不能用解析表达式表示的函数。数值逼近的基本方法是插值法。初等数学里 的三角函数表,对数表中的修正值,就是根据插值法制成的。

在遇到求微分和积分的时候, 如何利用简单的函数去近似代替所给的函数, 以便容易求微分和求积分, 也是计算方法的一个主要内容。微分方程的数值解法也是近似解法。常微分方程的数值解法有欧拉法、预 测校正法等。偏微分方程的初值问题或边值问题,目前常用的是有限差分法、有限元素法等。 有限差分法的基本思想是用离散的、只含有限个未知数的差分方程去代替连续变量的微分方程和定解 条件。求出差分方程的解法作为求偏微分方程的近似解。 有限元素法是近代才发展起来的,它是以变分原理和剖分差值作为基础的方法。在解决椭圆形方程边 值问题上得到了广泛的应用。目前,有许多人正在研究用有限元素法来解双曲形和抛物形的方程。 数值计算方法的内容十分丰富,它在科学技术中正发挥着越来越大的作用。

1.1.3 数值计算方法的特点
第一,面向计算机,要根据计算机特点提供实际可行的有效算法。即算法只能包括加、减、乘、除和 逻辑运算,是计算机能直接处理的。 第二,有可靠的理论分析,能任意逼近并达到精度要求,对近似算法要保证收敛性和数值稳定性,还 要对误差进行分析。 第三,要有好的计算复杂性,时间复杂性好是指节省时间,空间复杂性好是指节省存储量,这也是建 立算法要研究的问题,它关系到算法能否在计算机上实现。 第四,要有数值实验,即任何一个算法除了从理论上要满足上述三点外,还要通过数值实验证明是行 之有效的。 根据“数值计算方法”的特点,学习时我们要注意掌握方法的基本原理和思想,要注意方法处理的技巧 及其与计算机结合,要重视误差分析、收敛性及稳定性的基本理论;其次,要通过例子,学习使用各种数 值方法解决实际计算问题;最后,为了掌握本课的内容,还应做一定数量的理论分析与计算问题练习。由 于本课内容包括了微积分、代数、常微分方程的数值方法,以及高级程序设计的方法,读者必须掌握这几 门课的基本内容才能学好这一课程。

1.2

误差来源与误差分析

1.2.1 误差的来源
早在中学我们就接触过误差的概念,如在做热力学实验中,从温度计上读出的温度是 23.4℃就不是一 个精确的值,而是含有误差的近似值。事实上,误差的在我们的生活中无处不在,无处不有。如量体裁衣, 量与裁的结果都不是精确无误的,都有误差。人们可能会问:如果使用计算机来解决问题,结果还会有误 差吗?下面我们通过考察数学方法解决实际问题的主要过程来思考这个问题。 用数学方法解决一个具体的实际问题,首先要建立数学模型,这就要对实际问题进行抽象、简化,因 而数学模型本身总含有误差, 这种误差叫做模型误差。在数学模型中通常包含各种各样的参变量,如温度、 长度、电压等,这些参数往往都是通过观测得到的,因此也带来了误差,这种误差叫做观测误差。当数学 模型不能得到精确解时,通常要建立一套行之有效的数值方法求它的近似解,近似解与准确解之间的误差

就称为截断误差或方法误差。由于在计算机中浮点数只能表示实数的近似值,因此用计算机进行实际计算 时每一步都可能有误差,这种误差称为舍入误差。 例如,函数 f(x)用泰勒(Taylor)展开式

近似代替,则数值方法的截断误差是

( 介于 0 与 之间) 又如,在计算时用 3.14159 近似替代 -3.14159=0.0000026 产生的误差

就是舍入误差。

上述种种意味着都会影响计算结果的准确性,因此需要了解与研究误差。在数值计算方法中将着重研 究截断误差、舍入误差,并对它们的传播与积累作出分析。

1.2.2 绝对误差、相对误差与有效数字
1.绝对误差与绝对误差限

定义 1.1 设 为准确值,

为 x 的一个近似值,称

为近似值的绝对误差,或误差。

通常我们无法知道准确值 ,也不能算出误差的准确值 不超过某个正数 落在区间 ,则

,只能根据测量或计算估计出误差的绝对值 ,即

称为绝对误差限。有了绝对误差限就可知道 x 的范围 内。

例如用毫米测度尺测量一长度 ,读出的长度为 23mm,则有 对误差是有量纲和单位的。

mm。由此例也可以看到绝

2.相对误差与相对误差限
只用绝对误差还不能说明数的近似程度,例如甲打字时每百个字错一个,乙打字时平均每千个字错一 个,他们的误差都是错一个,但显然乙要准确些。这就启发我们除了要看绝对误差大小外,还必须顾及量 的本身。

定义 1.2 把近似值的误差

与准确值 的比值

称为近似值

的相对误差,记作



实际计算时,由于真值 通常是不知道的,通常取

。相对误差也可正可负,它的绝对值的上

界叫做相对误差限。记作

。即

。根据定义,甲打字时的相对误差

,乙打字时的

相对误差

。易知相对误差是一个无量纲量。

3.有效数字

当 有多位时,常常接四舍五入的原则得到 的前几位近似值

,例如

取3位



取5位



它们的误差都不超过末位数字的半个单位,即

现在我们将四舍五入抽象成数学语言,并引入一个新名词“有效数字”来描述它。

定义 1.3 若近似数 有 位有效数字。

的误差限是某一位的半个单位, 该位到

的第一位非零数字共有 位, 我们就说

如取 数字。



的近似值,

就有 3 位有效数字;取



的近似值,

就有 5 位有效

有效数字也可采用以下定义:

x 的近似数 x*写成标准形式:

(1.1)

其中

是 0 到 9 的一个数字,且

不为 0, 为整数,若

(1.2)

则称

有 位有效数字。

例 1.1 依四舍五入原则写出下列各数具有 5 位有效数字的近似数, 913.95872,39.1882,0.0143254,8.000033 : 解: 按定义,上述的 5 位有效数字的近似数分别为 913.96,39.188,0.014325,8.0000 注意,8.000033 的 5 位有效数字的近似值是 8.0000 而不是 8,8 只有一位有效数字。

4.有效数字与误差关系
不难看出,若由(1.1)给出某近似数有 位有效数字,则可以从(1.2)求得这个数的绝对误差限

因此,在 相同的情况下, 越大则

就越小,故有效数字越多,绝对误差限越小。

定理 1.1 用(1.1)表示的近似数

,若

具有 位有效数字,则其相对误差限为

反之,若

的相对误差限



至少具有 位有效数字。

证明 由式(1.1)得

所当

有 位有效数字时

反之,若

的相对误差限

有:

由式(1.2)式知道,

至少有 位有效数字,证毕

这说明有效数字越多,相对误差限越小。

例 1.2 要使

的近似值的相对误差限小于

,要取几个有效数字?

解:因为

,所以

,令

故取

即可满足。

5.函数计算的误差估计
一般情况,当自变量有误差时计算相应的函数值也会产生误差,其误差限可利用函数的泰勒展开式估 计。

设 开式

是一元函数, 的近似值为

,以

近似

,其误差界记作

,可用泰勒展

介于 ,

之间,取绝对值得

假定



的比值不太大,可忽略

的高阶项,于是可得计算函数值的误差限

例 1.3 设 解:

的相对误差限为 2%,则

的相对误差限是多少?

所以

的相对误差限为 0.02n。

当 近似值

为多元函数时,如计算 ,于是函数值 的误差

。如果

的近似值为

,则



由泰勒展开式得

于是误差限为

例 1.4 已测得某场地长 的值为

,宽

的值为

,已知

,试求面积的绝对误差限与相对误差限。

解:因

,那么

其中

于是绝对误差限

相对误差限

1.3

误差传播与若干防治办法

由前述可知,在数值计算中每步都可能产生误差。而一个问题的解决,往往要经过成千上万次运算, 我们不可能每步都有加以分析,下面,通过对误差的某些传播规律的简单分析,指出在数值计算中应注意 的几个原则,它有助于鉴别计算结果的可靠性并防止误差危害现象的产生。

1.要避免两相近数相减

在数值运算中两相近数相减会使有效数字严重损失,例如 字,但 只有两位有效数字,所以要尽量避免这类运算。

都具有五位有效数

通常采用的方法是改变计算公式,例如当



很接近时,由于

那么可用右端的

公式代替左端的公式计算,有效数字就不会损失。当 很大时,

也可用右端来代替左端。一般情况,当

时,可用泰勒展开

取右端的有限项近似左端。 如果计算公式不能改变,则可采用增加有效位数的方法。

“ ” 2.要防止大数“吃掉”小数
若参加运算的数的数量级相差很大,而计算机的位数有限,如不注意运算次序,就可能出现大数“吃 掉”小数的现象,影响计算结果。 例如在五位十进制计算机上,计算

写成规格化形式

由于计算时要对阶,

在计算机中表示为 0,计算出来

,结果严重失真!

如果计算时,先将 计算出来,再与 52492 相加,就不会出现大数“吃掉”小数的现象了。

3.注意简化计算步骤,减少运算次数
同样一个计算题,如果能减少运算次数,不但可节省计算机的计算时间,还能减少舍入误差。

例如,计算 x255 的值,如果逐个相乘要用 254 次乘法,但若写成 只要做 14 次乘法运算即可。

4.绝对值太小的数不宜作除数

设 与

分别有近似值

的近似值

,其绝对误差

显然,当

很小时,近似值

的绝对误差

有可能很大。因此,不宜把绝对值太小的数作除数。

第2章

非线性方程求根

与线性方程相比,非线性方程问题无论是从理论上还是从计算公式上,都要复杂得多。对于一般的非 线性方程 ,计算方程的根既无一定章程可寻也无直接法可言。例如,求解高次方程组 的根,求解含有指数和正弦函数的超越方程 或非线性方程组也是计算方法中的一个主题。一般地,我们用符号 般形式表示 为 ,方程的解称为方程的根或函数的零点。 的零点。解非线性方程 来表示方程左端的函数,方程一

通常,非线性方程的根不止一个,对于非线性方程,一般用迭代法求解。因此,在求解非线性方程时, 要给定初始值或求解范围。

2.1

实根的对分法

2.1.1 使用对分法的条件

对分法或称二分法是求方程近似解的一种简单直观的方法。设函数 , 则 在



上连续,且

上至少有一零点, 这是微积分中的介值定理, 也是使用对分法的前题条件。

计算中通过对分区间、缩小区间范围的步骤搜索零点的位置。

例 2.1 用对分法求解 解:

在区间

之间的根。

(1)



,由介值定理可得有根区间



(2)计算

,有根区间



(3)计算

(1.5+2)/2=1.75,

=0.078125,有根区间



(4)一直做到

(计算前给定的精度)或

时停止。

2.1.2 对分法求根算法

计算

的一般计算步骤如下:

1.输入求根区间

和误差控制量 ,定义函数



IF

THEN 退出选用其他求求根方法

2.WHILE



(1)计算中点 (2)分情况处理

以及

的值;

: 停止计算

,转向步骤 4

修正区间

修正区间

ENDWHILE

3.



4.输出近似根



图 2.1 给出对分法示意图

图 2.1 对分法示意图

在算法中, 常用

代替

的判断, 以避免

数值溢出。

对分法的算法简单,然而,若 使 在 上有零点,也未必有



是有几个零点时,只能算出其中一个零点;另一方面,即 。这就限制了对分法的使用范围。对分法只能计算方程

的实根。

2.2

迭代法

对给定的方程

,将它转达换成等价形式: ,如果迭代法收敛,即

。给定初值 ,有

,由此来构造迭代序列 ,则 为方程的根。 就是方程

的根。在计算中当

小于给定的精度控制量时,取

例如,代数方程

的三种等价形式及其迭代格式如下:

(1)

迭代格式

(2)

,迭代格式

(3)

迭代格式

对于方程 与迭代的初值有关?

构造的多种迭代格式

,怎样判断构造的迭代格式是否收敛?收敛是否

定理 2.1

若定义在

上,如果

满足

(1)当有



(2)



上可导,并且存在正数

,使对任意的

,有



则在 意的初值

上有惟一的点 均收敛于

满足 的不动点

,称



的不动点。而且迭代格式

对任

,并有误差估计式

(2.1)

证明:(1)先证明存在性:令

,则有

故有

,使得



(2)再证明惟一性:设

都是

的不动点,且

,则有

与假设矛盾,这表明

,即不动点是惟一的。

(3)当 定理

时,由于

可用归纳法证明,迭代序列 ,

,于是由微分中值



,得

(2.2)

因为

,所以当

时,

即迭代格式

收敛。

(4)误差估计式:

设 固定,对于任意的正整数

有,

由于

的任意性及

,故有

注:定理 2.1 是判断迭代法收敛的充分条件,而非必要条件。

要构造满足定理条件的等价形式一般是难于做到的。事实上,如果 形式 而 ,这时若初值 一,等价形式 函数 应满足 ,由 的连续性,一定存在 的邻域



的零点,若能构造等价 ,其上有

迭代也就收敛了。由此构造收敛迭代式有两个要素,其 ;其二,初值必须取自 。 的充分小邻域,这个邻域大小决定于

,及做出的等价形式

例 2.2 求代数方程

,在

附近的实根。

解:1)



,当

构造的迭代序列收敛。取



=2.08008,

=2.09235,

=2.094217

=2.094494,

=2.094543,

=2.094550

准确的解是 =2.09455148150 2)将迭代格式写为

,当

迭代格式

不能保证收敛。

2.3

迭代法的加工

对于收敛的迭代过程,只要迭代足够多次,就可以使结果达到任意的精度,但有时迭代过程收敛缓慢, 从而使计算量变得很大,因此迭代过程的加速是个重要的课题。 以下介绍一种埃特金(Aitken)方法。



方程,构造加速过程,算法如下:

(1)预测:

(2)校正:

(3)改进: 有些不收敛的迭代法,经过埃特金方法处理后,变得收敛了。

例 2.3 求方程



附近的根



解:若采用迭代公式: 法:

迭代法是发散的,我们现在以这种迭代公式为基础形成埃特金算



计算结果如表 2.1 所示: 表 2.1 计算结果

0 1 2 3 4 5

2.37500 1.84092 1.49140 1.34710 1.32518

13.3965 5.23888 2.31728 1.44435 1.32714

1.5 1.41629 1.35565 1.32895 1.32480 1.32472

我们看到,将发散的迭代公式通过埃特金方法处理后,竟获得了相当好的收敛性。

2.4

牛顿迭代法

对方程

可构造多种迭代格式

, 牛顿迭代法是借助于对函数

作泰勒展开

而构造的一种迭代格式。



在初始值

作泰勒展开:

取展开式的线性部分作为

的近似值,则有



,则



类似地,再将



作泰勒展开并取其线性部分得到:

一直做下去得到牛顿迭代格式:

,

(2.3)

牛顿迭代格式对应于

的等价方程是

(2.4)

若 是

的单根时, 的

则有

,只要初值

充分接近 ,有



所以牛顿迭代收敛。当 为

重根时,取下面迭代格式:

牛顿法的几何意义



为斜率作过

点的切线,即作



点的切线方程



,则得此切线与 轴的交点

,即

再作



处的切线,得交点

,逐步逼近方程的 根。如图 2.2 所示

图 2.2 牛顿切线法示意图

在区域 线性部分替代 。

的局部“以直代曲”是处理非线性问题的常用手法。在泰勒展开中,截取函数展开的

例 2.4 用牛顿迭代法求方程

,在

附近的根。

解: 计算结果列于表 2.2 中。 表 2.2 计算结果

K
0 1 2 3 4 5

xk
1.00 1.41176 1.62424 1.6923 1.69991 1.7

f (x )
-2.8 -0.727071 -0.145493 -0.0131682 -0.0001515 0

比较表 2.1 和表 2.2 的数值,可以看到牛顿迭代法的收敛速度明显快于对分法。

牛顿迭代法也有局限性。在牛顿迭代法中,选取适当迭代初始值

是求解的前题,当迭代的初始值

在某根的附近时迭代才能收敛到这个根,有时会发生从一个根附近跳向另一个根附近的情况,尤其在导数 数值很小时,如图 2.3 所示。

图 2.3 失效的牛顿迭代法

牛顿迭代算法

1.定义函数



,输入或定义迭代初始值

和控制精度 epsilon。

2.FOR

to MAXREPT

THEN{输出满足给定精度的近似解

,结束}

ENDFOR

3.输出:在

附近

无根。

注:在伪码中,“//”表示注释语句。

伪码又称过程设计语言,主要用于描述算法。它是某种高级语言和自然语言的混杂语言,它取某种高 级语言中的一些关键字,用于描述算法的结构化构造和数据说明等。伪码的语句中嵌有自然语言的叙述, 伪码易于

2.5
弦截法迭代格式

弦截法

在牛顿迭代格式中:

用差商 写成如下形式:

代替导数

,并给定两个初始值



,那么迭代格式可

(2.5) 上式称为弦截法。 用弦截法迭代求根, 每次只需计算一次函数值, 而用牛顿迭代法每次要计算一次函数值和一次导数值。 但弦截法收敛速度稍慢于牛顿迭代法。

弦截法的几何意义

做过两点 做过 和 ,如图 2.4 所示



的一条直线(弦),该直线与 轴的交点就是生成的迭代点 的一条直线, 是该直线与 轴的交点,继续做下去得到方程的根

, 再

图 2.4 弦截法示意图

例 2.5 用弦截法求方程

的根,取



解: 计算结果列于表 2.3 中。 表 2.3 计算结果

0 1 2 3 4 5 6 7

1.5 4 1.90909 1.65543 1.71748 1.70116 1.69997 1.7

-0.45 2.3 0.248835 -0.0805692 0.0287456 0.00195902 -0.0000539246 9.45910-8

弦截法算法

1.定义函数

,输入或定义控制精度值 epsilon 和迭代初始值



计算 2.FOR:=0,1…MAXREPT

表示

2.1

2.2

//

表示

2.3

THEN{输出满足给定精度的近似解 2.4

,结束}

//为下一次迭代准备数值 ENDFOR

3.输出:在初始值

附近

无根。

2.6

非线性方程组的牛顿方法

为了叙述简单,我们以解二阶非线性方程组为例演示解题的方法和步步骤,类似地可以得到解更高阶 非线性方程组的方法和步骤。 设二阶方程组

(2.6)

其中

为自变量。为了方便起见,将方程组写成向量形式:

,其中=



附近作二元泰勒展开,并取其线性部分,得到方程组



则有

如果

,解出

,得

再列出方程组:

解出

,得出

继续做下去,每一次迭代都是解一个方程组

(2.7)

即 例 2.6 求解非线性方程组

,直到

为止。

取初始值



解:

解方程得

继续做下去,直到

时停止。

将两个变量的非线性方程组推广到多变量的非线性方程组:

(2.8)



的向量形式: 用牛顿迭代法求方程组的解为:

其中:



为雅可比(jacobi)矩阵。计算中采用

(2.9)

一直做到

小于给定精度为止。



的领域中若

,则迭代收敛。

2.7

程序示例

程序 2.1 用牛顿迭代法求 算法描述



附近的根。

给定

,从

开始,根据牛顿迭代公式

计算



附近的根。

程序源码 ///////////////////////////////////////// // Purpose :牛顿迭代法求根 //

///////////////////////////////////////// #include<stdio.h #include<math.h>

//iterate(x)为 #define iterate(x) (x-((x*x-3)*x-1)/(3*x*x-3))

#define

((x*x-3)*x-1)

#define x0 #define

1.5

//迭代初值 x0

MAXREPT 1000 //最大迭代次数 //精度

#define epsilon 0.0001 void main( ) { int i;

double x_k=x0,x_k1=x0 for(i=i<MAXREPT;i++) { printf(“Got…%f/n”,x_k1); x_k1=iterate(x_k); //迭代

if fabs(x_k1 –x_k)<epsilon||fabs(f(x_k1)<epsilon { printf(“!Root:%f\n”,x_k1); //满足精度,输出 return; { x_k=x_k1 { printf(“After % d repeate, no solved.\m”,MAXREPT); //----------------------End of File----------------------计算实例

已知 程序输入输出

,取



,用牛顿迭代法计算

的根。

由本程序预设的 ! Root: 1.879385



,得到输出

程序 2.2 用弦截法求 算法描述



附近的根。

给定

,从

开始,根据弦截法迭代公式

求得 程序源码

在其附近的根。

///////////////////////////////////////// // Purpose:弦截法求根 //

///////////////////////////////////////// #include<stdio.h> #include<math.h> #define f(x) (x*x*x –7.7*x*x+19.2*x –15.3)//f 函数 #define x0 0.0 #define x1 1.0 #define MAXREPT 1000 #define epsilon 0.00001 void main( ) { int I; double x_k=x0,x_k1=x1,x_k2=x1; for(I=0;i<MAXREPT;i++= { printf(“!Root:%f\n”,x_k2); //最大迭代次数 //求解精度 //初值 x0,x1

x_k2=x_k1-(f(x_k1)*x_k1-x_k))/(f(x_k1)-f(x_k));//弦截求新 x_n if fabs(x_k2-x_k1) <epsilon||fabs(f(s_k2)) <epsilon { printf(“!Root:%f\n”,x_k2); //满足精度,输出 return; { x_k=x_k1; x_k1=x_k2; { printf(“After %d repdate, no solvde.\n”,MAXREPT); //-------------------End of File---------------------计算实例 //反复

用弦截法求方程 程序输入输出

的根,取

由本程序的 !Root: 1.700000



得到输出

第3章

解线性方程组的直接法


在近代数学数值计算和工程应用中,求解线性方程组是重要的课题。例如,样条插值中形成的

系式,曲线拟合形成的法方程等,都落实到解一个 元线性方程组,尤其是大型方程组的求解,即求线性 方程组(3.1)的未知量 的数值。

(3.1) 其中 ai j,bi 为常数。上式可写成矩阵形式 Ax = b,即

(3.2)

其中,

为系数矩阵,

为解向量, 的解存在且惟一,且有

为常数向量。当

detA=D

0 时,由线性代数中的克莱姆法则,方程组

为系数矩阵

的第 列元素以 代替的矩阵的行列式的值。 克莱姆法则在建立线性方程组解的理论

基础中功不可没,但是在实际计算中,我们难以承受它的计算量。例如,解一个 100 阶的线性方程组,乘 除法次数约为(101·100!·99),即使以每秒 的运算速度,也需要近 年的时间。在石油勘探、

天气预报等问题中常常出现成百上千阶的方程组,也就产生了各种形式方程组数值解法的需求。研究大型 方程组的解是目前计算数学中的一个重要方向和课题。 解方程组的方法可归纳为直接解法和迭代解法。从理论上来说,直接法经过有限次四则运算,假定每 一步运算过程中没有舍入误差,那么,最后得到方程组的解就是精确解。但是,这只是理想化的假定,在 计算过程中,完全杜绝舍入误差是不可能的,只能控制和约束由有限位算术运算带来的舍入误差的增长和 危害,这样直接法得到的解也不一定是绝对精确的。

迭代法是将方程组的解看作某种极限过程的向量极限的值,像第 2 章中非线性方程求解一样,计算极 限过程是用迭代过程完成的,只不过将迭代式中单变量 换成向量而已。在用迭代算法时,我

们不可能将极限过程算到底,只能将迭代进行有限多次,得到满足一定精度要求的方程组的近似解。 在数值计算历史上,直接解法和迭代解法交替生辉。一种解法的兴旺与计算机的硬件环境和问题规模 是密切相关的。一般说来,对同等规模的线性方程组,直接法对计算机的要求高于迭代法。对于中等规模 的线性方程组 ,由于直接法的准确性和可靠性高,一般都用直接法求解。对于高阶方程组和稀

疏方程组(非零元素较少),一般用迭代法求解。

3.1

消元法

3.1.1 三角形方程组的解
形如下面三种形式的线性方程组较容易求解。 对角形方程组

(3.3)



,对每一个方程,



显然,求解 n 阶对角方程的运算量为 下三角方程组



(3.4)

按照方程组的顺序,从第一个方程至第 个方程,逐个解出



由方程

,得

。将

的值代入到第二个方程





的值代入到第 个方程



计算

需要 次乘法或除法运算,

。因此,求解过程中的运算量为

上三角方程组

(3.5)

与计算下三角方程组的次序相反,从第 个方程至第一个方程,逐个解出



由第 个方程

。将

的值代入到第

个方程





的值代入到第 个方程

得解的通式

计算

需要

次乘法或除法运算

。因此求解过程中的运算量为

消元法的基本思想就是通过对方程组做初等变换,把一般形式的方程组化为等价的具有上述形式的易 解方程组。

3.1.2 高斯消元法与列主元消元法
高斯消元法
高斯消元法是我们熟悉的古老、简单而有效的解方程组的方法。下面是中学阶段解二元方程组(高斯 消元法)的步骤: (3.6) (3.7) 方程(3.6)乘以-3 加到第(3.7)个方程中得

代入(3.6)得



其方法相当于对方程组的增广矩阵做行的初等变换:

已是上三角矩阵,而

为原方程组的等价方程组,已化成易解的方程组形式。再用回代方法求解,得到: 这就是高斯消元法解方程组的消元和回代过程。 一般地,可对线性方程组(3.1)施行以下一系列变换; (1)对换某两个方程的次序; (2)对其中某个方程的两边同乘一个不为零的数; (3)把某一个方程两边同乘一个常数后加到另一个方程的两边。 记变换后的方程组为:

(3.8) 显然方程组(3.1)与(3.8)是等价方程组,或者说它们有相同的解。分别记方程组(3.1)与(3.8) 的增广矩阵为:

可以看出,

实际上是由

按一系列初等换后得到的

(1)对换

某两行元素;

(2)

中的某行乘一个不为零的数;

(3)把

的某一行乘一个常数后加到另一行。

高斯消元法就是通过以上(3)的变换,把

化为等价的上三角形式。

下面我们以 设方程组:

为例演示消元过程。

(3.9) 其增广矩阵为:

(1)若

,则将第一行乘以

加到第二行上;将第一乘以

加到第三行上;将第一行乘以

加到第四行上得到

(3.10)

即 其中:

(2)若 得到

则将第二行乘以

加到第三行上;将第二行乘以

加到第四行上,

(3.11) 其中:

(3)若

则将第三行乘以

加到第四行上,得到

(3.12)

其中:

已是上三角矩阵,于是得到了与原方程等价的易解形式的方程组:

(3.13)

再对方程组(3.13)依次回代解出



从式(3.12)可以得到系数矩阵

的行列式的值为

的对角元素的乘积。即

这也正是在计算机上计算方阵

的行列式的常规方法。

要将上述解方程步骤推广到 阶方程组,只需将对控制量“4”的操作改成对控制量 的操作即可。 元方程组高斯消元法的步骤如下:

对于

约定



(3.14)

经过以上消元,我们已得到与 见,仍记 的元素为

等价的方程组

,其中

已是一个上三角矩阵。为简单起

(3.15) 即已得到原方程组的解。

高斯消元法算法

在算法中略去了变量,矩阵和向量的定义部分。在消元过程中,将 1.输入:方程组阶数 n,方程组系数矩阵 A 和常数向量 b。 2.FOR k:=1 TO n-1 { FOR i:=k+1 TO n //消元过程

仍放在

元素的位置上。

{ FOR j:=k+1 TO n

// 假定

{

}

//

ENDFORJ

} }

// //

ENDFORI ENDFORK

3.FOR i:=n TO1 { s:=bi

//回代求解

FOR j:=i+1 TO n DO

}

4.输出方程组的解



高斯消元法的运算量
整个消元过程即式(3.14)的乘法和除法的运算量为

回代过程即式(3.15)的乘除运算量为

故高斯消元法的运算量为

(3.16)

高斯消元法的可行性
在上面的消元法中,未知量是按照在方程组中的自然顺序消去的,也叫顺序消元法。

在消元过程中假定对然元素 故高斯消元法可行的充分必要条件为

,消元步骤才能顺利进行,由于顺序消元不改变 的各阶主子式不为零。但是,实际上只要

的主子式值, 方程组

就有解。故高斯消元法本身具有局限性。

另一方面,即使高斯消元法可行,如果

很小,在运算中用它作为除法的分母

,会导致其它元素数量级的严重增长和舍入误差的扩散。这是高斯消元法的另 一缺陷。 (3.17)

例 3.1 方程组

(3.18)

的精确解为:

。在高斯消元法计算中取 5 位有效数字。

解:方程(3.17)×(-1)/0.0003+方程(3.18)得:

,代入方程(3.17)得 得到等价方程组

。由此得到的解完全失真,如果交换两个方程的顺序,

经高斯消元后有

由此可看到,在有些情况下,调换方程组的次序对方程组的解是有影响的,对消元法中抑制舍入误差 的增长是有效的。 如果不调换方程组的次序,取 6 位有效数字计算方程组的解,得到

取 9 位有效数字计算方程组的解,得到

由此可见有效数字在数值计算中的作用。

列主元消元法
如果在一列中选取按模最大的元素,将其调到主干方程位置再做消元,则称为列主元消元法。调换方 程组的次序是为了使运算中做分母量的绝对值尽量地大,减少舍入误差的影响。 用列主元方法可以克服高斯消元法的额外限制,只要方程组有解,列主元消元法就能畅通无阻地顺利 求解,同时又提高了解的精确度。

更具体地,第一步在第一列元素中选出绝对值最大的元素 的所有元素,再做化简 为零的操作。

,交换第一行和第



对于每个 对 行和

在做消元前,选出

中绝对值最大的元素 ,可证



行交换后,再做消元操作,这就是列主元消元法的操作步骤。由于

中至少有一个元素不为零,从理论上保证了列主元消元法的可行性。列主元消元法 与高斯消元法相比,只增加了选列主元和交换两个方程组(即两行元素)的过程。

列主元消元法算法

1.输入:方程组阶数 ,方程组系数矩阵

和常数向量项 。

2.

//选主元的消元过程

{//选择

// 交换第 行和第



} } 3.FOR i:=n TO 1

// ENDFORI // ENDFORK // 回代求解

4.输出方程组的解



如果对于第 步,从 行至 行和从 列至 列中选取按模最大的

,对第 行和第 行交换,对第 列和第 v 列交换, 这 就是全主元消元法。在 k 列和第 v 列交换时,还要记录下 v 的序号,以便恢复未知量 xk 和 xv 的位置。

3.1.3 高斯-若尔当(Gauss-Jordan)消元法
高斯消元法将系数矩阵化为上三角矩阵,再进行回代求解;高斯-若尔当消元法是将系数矩阵化为对 角矩阵,再进行求解,即对高斯消元法或列主元消元法得到的等价增广矩阵:

用第 4 行乘 到第 1 行上,得到

加到第 3 行上, 用第 4 行乘

加到第 2 行上, 用第 4 行乘



用第 3 行乘 加到第 1 行上,得到

加到第 2 行上, 用第 3 行乘

加到第 1 行上, 再用第 2 行乘

(3.19)

为方便起见,仍记式(3.19)系数矩阵元素为

,常数项元素为

。那么

用初等变换化系数矩阵为对角矩阵的方法称为高斯-若尔当消元法。从形式上看对角矩阵(3.15)比 上三角矩阵(3.11)更为简单,易于计算 略大于上三角矩阵回代的工作量。 ,但是将上三角矩阵(3.11)化为对角矩阵(3.15 )的工作量

— 高斯—若尔当消元法求逆矩阵



为非奇异矩阵,方程组 ,则 。

的增广矩阵为

。如果对 应用高斯-若尔当消元法化为

例 3.2 用高斯-若尔当消元法求

的逆矩阵



解:

解得:

3.2

直接三角分解法

仍以

为例,在高斯消元法中,从对方程组进行初等变换的角度观察方程组系数矩阵的演变过程。 和 。记

第一次消元步骤将方程组(3.9)化为方程组(3.10),相当于用了三个初等矩阵左乘 ,容易验证,



得到

其中

将方程组(3.10)化为方程组(3.11),相当于用了两个初等矩阵左乘

和 。





·

同理,将方程组(3.11)化为方程组(3.12),相当于用一个初等矩阵左乘

和 。



,有

完成了消元过程,有

亦有所有消元步骤表示为: 仍为下三角矩阵,并有

左乘一系列下三角初等矩阵。容易验证,这些下三角矩阵的乘积

于是有



这里

仍为下三角矩阵,其对角元素为 1,称为单位下三角矩阵,而

已是上三角矩阵。



,则有

结果表明若消元过程可行,可以将 解方程组的直接分解法。

分解为单位下三角矩阵

与上三角矩阵

的乘积。由此派生出

由高斯消元法得到启发,对 程。如果直接分解 出 ;再由 得到 和

消元的过程相当于将 , 。这时方程

分解为一个上三角矩阵和一个下三角矩阵的过 化为 ,令 ,由 解

,解出 。这就是直接分解法。

将方阵 分解;当

分解为

,当

是单位下三角矩阵,

是上三角矩阵时,称为多利特尔(Doclittle) 的

是下三角矩阵,

是单位上三角矩阵时,称为库朗(Courant)分解。分解的条件是若方阵

各阶主子式不为零,则多利特尔分解或库朗分解是可行和惟一的。

3.2.1 多利特尔分解
多利特尔分解步骤



的各阶主子式不为零,

可分解为

,其中

为单位下三角矩阵,

为上三角矩阵,即

(3.20)

矩阵



共有

个未知元素,按照

的行元素 或

的列元素的顺序,对每个 和

列出式(3.16)两边

对应的元素关系式,一个关系式解出一个

的元素。下面是计算

的元素的步骤:

(1)计算

的第一行元素

要计算

,则列出式(3.20)等号两边的第 1 行第 1 列元素的关系式:



。一般地,由

的第一行元素的关系式

得到

(2)计算

的第一列元素

要计算

,则列出式(3.20)等号两边的第 2 行第 1 列元素的关系式:



。一般地,由

的第 1 列元素的关系式

得到

(3)计算

的第 2 行元素

(4)计算

的第 2 列元素

……

若已算出

的前

行,

的前

列的元素,则

(5)计算

的第 行元素



的第 行元素的关系式:

是上三角矩阵,

列标

行标 。

得到

(3.21)

(6)计算

的第 列元素



的第 列元素的关系式:



是下三角矩阵,∴行标标

行标 。

得到

(3.22)

一直做到

的第

列,

的第 行为止。

用 同。

直接分解方法求解方程组所需要的计算量仍为

, 和用高斯消元法的计算量基本相

可以看到在分解中 间,可将 以及

的每个元素只在式(3.21)或(3.22)中做而且仅做一次贡献,如果需要节省空 相应元素的位置上。

的元素直接放在矩阵

用直接分解法解方程 三角矩阵,容易得到

,首先作出分解



,解方程组

。由于

是单位下

(3.23)

再解方程组

,其中

(3.24)



进行

分解时,并不涉及常数项 。因此,若需要解具有相同系数的一系列线性方程组时,使

用直接分解法可以达到事半功倍的效果。 多利特尔直接分解算法

1.输入:方程组阶数 ,系数矩阵

和常数项 。

2.

// 计算

的第 行元素

//计算

的第 列元素

}

//

完成

分解

3.

// 解方程组

4.

// 解方程组

5.输出方程组的解 例 3.2 用多利特尔分解求解方程组:

解:



,有



,有



,有



,即



,即

3.2.2 追赶法
很多科学计算问题中,常常所要求解的方程组为三对角方程组

(3.25) 其中

(3.26) 并且满足条件:

称 值函数的

为对角占优的三对角线矩阵。其特征是非零元素仅分布在对角线及对角线两侧的位置。在样条插 关系式中就出现过这类矩阵。事实上,许多连续问题经离散化后得到的线性方程组,其系数

矩阵就是三对角或五对角形式的对角矩阵。

利用矩阵直接分解法,求解具有三对角线系数矩阵的线性方程组十分简单而有效。现将 角矩阵 ,及单位上三角矩阵 的乘积。即 ,其中

分解为下三



(3.27)

利用矩阵乘法公式,比较

两边元素,可得到

于是有

(3.28)

由些可见将

分解为



,只需计算



两组数,然后解

,计算公式为:

(3.29)

再解

则得

(3.30)

整个求解过程是先由(3.28)及(3.29)求 求出 ,这时





, 这时

是“追”的过程, 再由(3.24)

是往回赶,故求解方程组(3.25)的整个过程称为追赶法。

3.2.3 对称矩阵的 LDLT 分解
对称正定矩阵也是很多物理问题产生的一类矩阵,正定矩阵的各阶主子式大于零。可以证明,若 定,则存在下三角矩阵 ,使 ,直接分解 正

的分解方法,称为平方根法。由于在平方根 分

法中含有计算平方根的步骤,因此很少采用平方根法的分解形式。对于对称矩阵,常用的是 解。



作多利特尔分解

,即

(提出

矩阵的对角元素)



对称正定,可得

,令



的对称性和分解的惟一性可证





(3.31)

是对角元素为 1 的单位下三角矩阵。

对矩阵

作多利特尔或库朗分解,共计算

个矩阵元素;对称矩阵的

分解,只需计算 分解计算公式。

个元素,减少了近一半的工作量。借助于多利特尔或库朗分解计算公式,容易得到



有多利特尔分解形式:

其中

在分解可套用多利特尔分解公式,只要计算下三角矩阵 的元素, 的 行 列的元素用 的



的对角元素

。计算中只需保存

表示。由于对称正定矩阵的各阶主子式大于零,直接调用

多利特尔或库朗分解公式可完

分解计算,而不必借助于列主元的分解算法。

对于



(3.32)

(3.33)



,解方程组

可分三步完成:

(1)解方程组

,其中



(3.34)

(2)解方程组

其中



(3.35)

(3)解方程

(3.36)

对称矩阵的

分解算法

1.输入:方程组阶数

,系数矩阵

和常数项 。

2.

3.略去解方程组步骤 从矩阵分解角度看,直接分解法与消元法本质上没有多大区别,但实际计算时它们各有所长。一般来 说,如果仅用单字长进行计算,列主元消元法具有运算量较少、精度高的优点,故是常用的方法。但是, 为了提高精度往往采取单字长数双倍内积的办法(即作向量内积计算时,采用双倍加法,最终结果再舍入 成单字长数),这时直接分解法能获得较高的精度。

例 3.4 用

分解求解方程组:

解:



,有

3.3





3.3.1 向量范数

在一维空间中,实轴上任意两点

距离用两点差的绝对值

表示。绝对值是一种度量形式的定义。

范数是对函数、向量和矩阵定义的一种度量形式。任何对象的范数值都是一个非负实数。使用范数可 以测量两个函数、向量或矩阵之间的距离。向量范数是度量向量长度的一种定义形式。范数有多种定义形 式,只要满足下面的三个条件即可定义为一个范数。同一向量,采用不同的范数定义,可得到不同的范数 值。

定义 3.1 对任一向量 足下面三个性质: (1) ,有

,按照一个规则确定一个实数与它对应,记该实数记为

,若



,当且仅当

时,

(非

(3.37)

负性) (2) (3) , , ,有 ,有 (齐次性) (三角不等式)

那么称该实数

为向量

的范数。

几个常用向量范数

向量



范数定义为

其中,经常使用的是

三种向量范数。

或写成

例 3.5 计算向量

的三种范数。

向量范数的等价性

有限维线性空间 义,则必存在

中任意向量范数的定义都是等价的。若 ,使 均有



上两种不同的范数定

或 (证明略)

向量的极限
有了向量范数的定义 ,也就有了度量向量距离的标准,即可定义向量的极限和收敛概念了。



为 是收敛的(

上向量序列,若存在向量 是某种向量范数),

使

,则称向量列

称为该向量序列的极限。

由向量范数的等价知,向量序列是否收敛与选取哪种范数无关。

向量序列 即 存在。



收敛的充分必要条件为其序列的每个分量收敛,



,则

就是向量序列

的极限。

例 3.6 求向量序列 解:算出每个向量分量的极限后得

极限向量。

在计算方法中,计算的向量序列都是数据序列,当 向量 。

小于给定精度时,取

为极限

3.3.2 矩阵范数
矩阵范数定义

定义 3.2 如果矩阵

的某个非负实函数

,记作

,满足条件:

(1)

当且仅当

时,

(非负性)

(2)

(齐次性)

(3)对于任意两个阶数相同的矩阵



(三角不等式性)

(4)

矩阵为同阶矩阵

(相容性)

则称

为矩阵范数。

矩阵的算子范数
常用的矩阵范数是矩阵的算子范数,可用向量范数定义:



,记方阵

的范数为

,那么



(3.38)

称为矩阵的算子范数或从属范数。这样定义的矩阵范数满足矩阵范数的所有性质外,还满足相容性:

为 阶矩阵,

恒有

(3.39)

根据定义,对任一种从属范数有

,即单位矩阵的范数是 1。

常用矩阵范数
向量有三种常用范数,相对应的矩阵范数的三种形式为:



的行范数)

(3.40)



的列范数)

(3.41)





的最大特征值)(

的 2 范数)

(3.42)

证明:既然矩阵的算子范数是 就找到了矩阵的范数。

上满足

向量范数

的上确界,那么,找到这个上确界也

(1)任取

,则



另一方面设极大值在

列达到,即



, 除第 个

分量为 1 外,其余分量均为 0,于是有



由于

,故

因此有

(2)任取 即



,则

另一方面设极大值在 行达到,取 这里

于是



(3)

为对称非负矩阵,具有非负特征值,并具有 个相互正交的单位特征向量。



的特征值为

,相应的特征向量为

,其中

为相互正交的单位向量。 设

,并且

,即

,则

即对任意

均有

, 故



,则有



于是

如果 A 是对称矩阵,那么

,设

的特征值是 , 则有

还有一种与向量范数 表示,其定义为

相容的矩阵范数, 称为欧几里得 (Euclid) 范数或舒尔 (Schur) 范数, 用

(3.43)

因为欧几里得范数易于计算,在实用中是一种十分有用的范数。但它不能从属于任何一种范数,因为 。 与向量范数的等价性质类似,矩阵范数之间也是等价的。

例 3.7 : 解:

,分别有



的特征值

矩阵范数等价性

定理 对

上的任两种范数



,存在常数



,使



t≧

矩阵特征值与范数关系



是矩阵

的特征值(即存在非零向量

使得:

),对任一算子范数有



(相容性)

即矩阵特征值的模不大于矩阵的任一范数。

谱半径



有特征值

记 。

则称



的谱半径。有了谱半径的定义,矩阵的 2

范数可记为:

谱半径与矩阵范数关系

由矩阵谱半径定义,可得到矩阵范数的另一重要性质,



3.4

矩阵的条件数

在解方程组时,我们总是假定系数矩阵

和常数项 是准确的,而在实际问题中,系数矩阵

和常数

项 往往是从前面的近似计算所得,元素的误差是不可避免的。这些误差会对方程组 大的影响?矩阵的条件数将给出一种粗略的衡量尺度。

的解 产生多

定义 3.3 若

非奇异称



的条件数。其中

表示矩阵的某种范数。

用矩阵

及其逆矩阵

的范数的乘积表示矩阵的条件数,由于矩阵范数的定义不同,因而其条件数

也不同,但是由于矩阵范数的等价性,故在不同范数下的条件数也是等价的。矩阵条件数的大小是衡量矩 阵“坏”或“好”的标志。

对于线性方程组

,若系数矩阵有小扰动

这时方程组的解也有扰动,于是





的影响可表示为:

(3.44)

若常数 有小扰动

,其对

的影响表示为:

(3.45)

因此,

大的矩阵称为“坏矩阵”或“病态矩阵”。对于

大的矩阵,小的误差可能会引

起解的失真。一般说来,若 地,当 很小时,

的按模最大特征值与按模最小特征值之比值较大时,矩阵就会呈病态。特别

总是病态的。

例 3.8 方程组

方程组的解为

对常数项 引入扰动

,则解为



。虽然

很小,但解已完全失真。



的病态是显而易见的。

3.5

病态方程组的解法

如果系数矩阵 下措施:

的条件数

远大于 1,则

为病态方程组,对病态方程组求解可采用以

(1)采用高精度运算,减轻病态影响,例如用双倍字长运算;

(2)用预处理方法改善 价,而 时可选择 的条件数比

的条件数,即选择非奇异阵 改善,而求 的解

使 即





为原方程组的解,计算

为对角阵或三角阵;

(3)平衡方法,当

中元素的数量级相差很大,可采用行平衡或列平衡的方法改善

的条件数。设

非奇异,计算 价于求 ,或 这时



,于是求



的条件数可得到改善,这就是行均衡法。

例 3.9 给定方程组



A 的条件数

,若用行均衡法可取

,则平衡后的方程



,用三位有效数字的列主元消元法求解得



3.6
程序 3.1 用列主元高斯消元法求解方程组:

程序示例

算法描述
输入:方程组阶数 n,矩阵 a 及列向量 b。

分解过程:对

{

选择

交换第 行和第

行方程





回代过程:对

解得

程序源码
//////////////////////////////////////

// ////////////////////////////////////// #include<stdio.h> #include<stdio.h> #defineMAX-N20 int main() { int n; int i,j,k; int mi;double mx, tmp; static ouble a[MAXN] ,[MAXN],b[MAXN],x[MAX N]; printf(“\ninput nvalue (dim of Ax=b):”); //输入 Ax=b 的维数 scanf(“%d”,&n); if(n>MAXN) { printf(“The input n is larger than MAXN,please redefine the MAXN.\n”); return 1; } if(n<=0) { printf(“Please input a number between 1 and % d.\n”,MAX_N); return 1; { //输入 Ax=b 的 A 矩阵 //(xi,yi)的最大维数

//

printf(“Now input the matrix a(i, j), i, j=0,…,% d: \n”, n-1); for(i=0, i<n; j++) for(j=0;j<n;j++ scanf(“%1f”,&a[i][j]) //输入 b 向量 printf(“ Now input the matrix b(i), i=0,…,% d: \n”, n-1); for(i=0;i<n; i++)scanf(“% 1f ”, &b[i]); for(i=0;i<n-1; i++) { for(j=i+1,mi=i, mx=fabs(a[i][i]); j<n;j++) //找主元素 if(fabs(a[j][i])>mx) { mi=j; mx=fabs(a[j][i]); { if(i<mi) { tmp=b[i]; b[i]=b[mi]; b[mi]=tmp; for(j=I; j<n; j++) { tmp=a[i][j]; a[i][j]=a[mi][j]; a[mi][j]= tmp; } } //交换两行

for(j=i+1; j<n; j++) { tmp=-a[j][i]/a[i][i]; b[j]+=b[i]*tmp; for=(k=i, k<n; k++) a[j][k]+=a[i][k]*tmp; }} x[n-1]=b[n-1]/a[n-1][n-1]; for(i=n-2; i>0; i--) { x[i]=b[i]; for(j=i+1, j<n; j++) x[i]- =a[i][j]*x[j]; x[i]/=a[i][i]; } printf(“Slove…x_i=\n”);

//高斯消元

//求解方程

//输出

for(i=0; i<n++)printf(“%f\n,x[i]”); return 0; } //……………………End of File……………………

计算实例
用列主元高斯消元法水求解方程组:

程序输入输出
Input n value(dim of Ax=b): 3

Now input the matrix a(i, j), i, j=0, …, 2: 2 1 1 1 3 2 1 2 2 Now input the matrix b(i), i=0, …, 2: 4 6 5 Solve… x_i= 1.000000 1.000000 1.000000 程序 3.2 用库朗分解公式求解方程组:

算法描述

输入方程组阶数 , 矩阵

及列向量 ;

将矩阵

分解为

,其中











程序源码
//////////////////////////////// // Purpose:库朗分解解 Ax=b // /////////////////////////////// #include<stdio.h> #include<stdio.h> #define MAXN20 int main( ) { int n; int i, j, k; int mi;double mx, tmp; statc doudle a[MAXN] [MAXN],b[MAXN],x[MAX N],y[MAXN]; statc double I[MAXN] [MAXN],u[MAXN] [MAXN]; //(xi,y i)的最大维数

printf(“nInput n value(dim of Ax=b):”);//输入

的维数

scanf(“%d”,&n); if(n>MAXN) { printf(“The input n is larger then MAXN,pleadefine the MAXN.\n”); return 1; if(n<=0) { printf(“Please input a number between 1 and %d.\n”, MAX N); retun 1; }

//输入



矩阵

printf(“Now input the matrix a(i,j),i,j=0…,%d.\n”,n-1); for(i=0;i<n;i++) for(j=0;j<n;j++) scanf(“%If”,&a[i][j]); //输入 b 矩阵 printf(“Now input the matrix b(i),i=0,…d:\n”,n-1); for(i=0;i<n;i++)scanf(“%If”,&b[i]); for(i=0;i<n;i++)u[i][i]=1; for(k=0;k<n;k++) //U 矩阵对角元素为 1

{ for(i=k;i<n;i++)

//计算

矩阵的第

列元素

{ 1[i][k]=a[i][k]; For(j=0;j<=k-1;j++, 1[i][k] -=(1[i][j]*u[j][k]); }

for(j=k+1;j<n; j++) { u[k][j]=a[k][j]; for(i=0; i<=k-1; i++)

//计算

矩阵的第

行元素

u[k][j] - =(1[k][j]*u[i][j]); u[k][j]/=1 [k][k]; } } for(i=0; i<n; i++) { y[i]=b[i]; for(j=0; i<=i–1; j++) y[i]-=(1[i][j]*y[j]); y[i]/=1[i][i]; { for(i=n-1; i>=0; i--) { x[i]=y[i]; //Ux=y 计算 x //Ly=b 计算 y

for(j=i+1; j<n; j++) x[i]-=(u[i][j]*x[j]); } printf(“Solve…x_I=\n”); //输出

for(i=0; i<=n; I++)printf(“% f\n”,x[i]); return 0; }

计算实例
用库朗分解公式求解方程组:

程序输入输出
Input n value(dim of Ax=b):3 Now input the matrix a(i,j)i,j=0,…,2: 1 2 1 -2 -1 -5 0 -1 6 Now input the matrix b(i),i=0,…,2: 24 –63 50 solve…x-i= 7.000000 4.000000 9.000000

程序 3.3 用

分解求解方程组:

算法描述

输入矩阵

及列向量 。

将矩阵

分解为

,其中





,先解

再由



最后由





程序源码
///////////////////////////////

// Purpose:

分解解方程组//

/////////////////////////////// #include<stdio.h> #include<math.h> #define MAX-N20 int main( ) { int n; //(x-i,y-i)的最大维数

int i, j, k; iInt mi;double mx,tmp; static double a[MAX-N][MAX-N],B [MAX-N],x[MAX-N],y[MAX-N],z[MAX-N]; static double L[MAX-N] [MAX-N],d[MAX-N]; printf(“\ninput n value(dim of Ax=b):”);//输入 Ax=b 的维数 scanf(“%d”,&n); If(n>MAX-N) { pritf(“The input n is lager then MAX-N,pleafine redefine the MAX-N.\n”); return 1; } if(n<=0) { printf(“Please input a numbr between 1 and %d.\n”, MAX-N); retun 1 } //输入 Ax=b 的 A 矩阵 printf(“Now input the matrix a(i, j), i, j=0, …, % d: \n”, n-1); for(i=0; i<n; i++)for (j=0;j<n;j+=scanf(“%1f ”, &b[i]); for(i=0; i<n; i++)L[i][i]=1; //L 矩阵对角元素为 1 for(k=0; k<n; k++) { d[k]=a[k][k]; for(j=0; j<=k-1; j++) d(k)-=(L[k][j]*d[j]); 计算 d_i

for(i=k+1;i<n; i++)

//计算

矩阵的第

列元素

{ L[i][k]=a[i][k]; For(j=0;j<=k-1;j++) L[i][k]-=(L[i][j]* L[i][j]*d[j]); L[i][k]/=d[k]; } }

for(i=0;i<n;i++) { z[i]=b[i]; for(j=0;j<=i-1;j++) z[i]-= L[i][j]*z[j]; }

//求解

的z

for(i=0;i<n;i++)

//求解



y[i]=z[i]/d[i];

for(i=n-1;i>=0 ; i--) { x[i]=y[i]; for(j=i+1;j<n;j++) x[i]-=[j][i]*x[j]; } printf(“Solve…x-i=\n”)

//求解

//输出

for(i=0;i<n;i++=printf(“%f\n”,x[i]);

return 0;

计算实例



分解求解方程组:

程序输入输出
Input n value(dim of Ax=6):3 Now input the matrix a(i,j),i,j,=0,…,2: 1 –11 –13 –2 1 –2 4.5

Now input the matrix b(i),i=0,…,2: 4 –8 12 Solve…x-i= 1.000000 -1.000000 2.000000

第4章
用迭代法求解线性方程组 变换,构造同解方程组 则由 得到

解线性方程组的迭代法
与第 4 章非线性方程求根的方法相似,对方程组 (对 可构造各种等价方程组,如分解 ),以此构造迭代关系式 进行等价 , 可逆,

(4.1)

任取初始向量

,代入迭代式中,经计算得到迭代序列



若迭代序列

收敛,设

的极限为

,对迭代式两边取极限



是方程组

的解,此时称迭代法收敛,否则称迭代法发散。我们将看到,

不同于非线性方程的迭代方法,解线性方程组的迭代收敛与否完全决定于迭代矩阵的性质,与迭代初始值 的选取无关。迭代法的优点是占有存储空间少,程序实现简单,尤其适用于大型稀疏矩阵;不尽人意之处 是要面对判断迭代是否收敛和收敛速度的问题。

可以证明迭代矩阵的与谱半径 的特征根。事实上,若 为方程组

是迭代收敛的充分必要条件,其中 的解,则有

是矩阵

再由迭代式

可得到

由线性代数定理,

的充分必要条件



因此对迭代法(4.1)的收敛性有以下两个定理成立。

定理 4.1 迭代法

收敛的充要条件是



定理 4.2 迭代法

收敛的充要条件是迭代矩阵

的谱半径

因此,称谱半径小于 1 的矩阵为收敛矩阵。计算矩阵的谱半径,需要求解矩阵的特征值才能得到,通 常这是较为繁重的工作。但是可以通过计算矩阵的范数等方法简化判断收敛的工作。前面已经提到过,若 ||A|| p 矩阵 和 的范数,则总有 。因此,若 ,则 必为收敛矩阵。计算矩阵的 1 范数

范数的方法比较简单,其中

于是,只要迭代矩阵 的。

满足



,就可以判断迭代序列

是收敛

要注意的是,当



时,可以有

,因此不能判断迭代序列发散。

在计算中当相邻两次的向量误差的某种范数 为方程组

小于给定精度时,则停止迭代计算,视

的近似解(有关范数的详细定义请看 3.3 节。)

4.1

雅可比(Jacobi)迭代法

4.1.1 雅可比迭代格式
雅可比迭代计算
元线性方程组

(4.2)

写成矩阵形式为

。若

将式( 4.2)中每个方程的 则得到下列同解方程组:

留在方程左边,其

余各项移到方程右边;方程两边除以



,构造迭代形式



(4.3)

迭代计算式(4.3)称为简单迭代或雅可比迭代。任取初始向量 列

,由式(4.3)可得到迭代向量序

雅可比迭代矩阵





,得到等价方程:



不难看出,

正是迭代式(4.3)的迭代矩阵,

是常数项向量。于是式(4.3)可写成矩阵形式:

(4.4) 其中:

雅可比迭代算法

下面描述解线性方程组 迭代要求,即

的雅可比迭代算法,为了简单起见,在算法中假定矩阵 ,并设由系数矩阵 构造迭代矩阵 是收敛的。

满足雅可比

1.定义和输入系数矩阵

与常数项向量

的元素。

2.FOR i:=1,2,…,n

{

//假定 FOR j:=1,2,…,n

,形成常数项向量

}

//形成迭代矩阵元素

3.

// 赋初始值,x1 和 x2 分别表示



4.WHILE

x1:=x2 x2:=B*x1+g // FOR u:=1,2,…,n
//

s:= g[u]; s:=s+b[u][v]*x1[v];

// FOR v:=1,2,…,n // ENDWHILE

x2[u]:=s;

5.输出方程组的解 例 4.1 用雅可比方法解下列方程组:

解:方程的迭代格式:



雅可比迭代收敛。

取初始值

,计算结果由表 4.1 所示。 表 4.1 计算结果

0 1 2 3 4 5 6 7

1 -1.5 -1.25 -0.915 -0.9575 -1.01445 -1.00722 -0.997543

1 1.6 2.08 2.068 1.9864 1.98844 2.00231 2.00197

1 0.9 1.09 1.017 0.9847 0.99711 1.0026 1.00049 0.25 0.48 0.335 0.0816 0.05695 0.01387 0.009687

方程组的准确解是

4.1.2 雅可比迭代收敛条件

对于方程组

,构造雅可比迭代格式

其中



。当迭

代矩阵的谱半径 时, 迭代收敛, 这是收敛的充分必要条件。 迭代矩阵的某范数 时, 迭代收敛。要注意的是范数小于 1 只是判断迭代矩阵收敛的充分条件,当迭代矩阵的一种范数||B||>1,并

不能确定迭代矩阵是收敛还是发散。 例如, 和 0.8。 是收敛矩阵。

, 则

, 但它的特征值是 0.9

当方程组的系数矩阵

具有某些性质时,可直接判定由它生成的雅可比迭代矩阵是收敛的。

定理 4.3 若方程组

的系数矩阵

,满足下列条件之一,则其雅可比迭代法是收敛的。

(1)

为行对角占优阵,即

(2)

为列对角占优阵,即

证明:(1)雅可比迭代矩阵

其中

(2)

为列对角优阵,故

为行对角占优阵,由系数矩阵

构造的迭代矩阵

为行

对角占优阵,则有



得到







由系数矩阵

构造的雅可比迭代矩阵收敛。

(如矩阵

既是行对角占优阵,也是列对角占优阵)

定理 4.4 若方程组系数矩阵

为对称正定阵,并且

也为对称正定,则雅可比迭代收敛。

4.2
高斯-赛德尔迭代计算

高斯-塞德尔(Gauss-Seidel)迭代法

在雅可比迭代中, 用 的计算公式是

的值代入方程 (4.2) 中计算出

的值,

事实上,在计算

前,已经得到

的值,不妨将已算出的分量直接代入迭代式中, 及

时使用最新计算出的分量值。因此

的计算公式可改为:

即用向量 用向量

计算出 计算出

的值,用向量

计算出

的值



的值,这种迭代格式称为高斯—塞德尔迭代。

对于方程组 AX=y ,如果由它构造高斯-塞德尔迭代和雅可比迭代都收敛,那么,多数情况下高斯— 塞德尔迭代比雅可比迭代的收敛效果要好,但是情况并非总是如此。

构造方程组 中每个方程的 组:

的高斯-塞德尔迭代格式的步骤与雅可比类似,设 留在方程的左边,其余各项都移到方程的右边;方程两边除以

将式(4.1) ,得到下列同解方程



,对方程组对角线以上的

取第 步迭代的数值,对角线以下的

取第

步迭代的数值,构造高斯—塞德尔迭代形式:

(4.5) 例 4.2 用高斯-塞德尔方法解方程组:

解:方程的迭代格式:

取初始值



时,

时,

计算结果如表 4.2 所示。 表 4.2 计算结果

0 1 2 3 4
— 高斯—塞德尔迭代矩阵

0 -2.5 -0.88 -1.0042 -1.0005

0 2.1 2.004 1.9984 2.0002

0 1.14 0.9876 1.0006 1.0000 2.5 1.62 0.1242 0.0037



写成等价矩阵表达式:

构造迭代形式:

有 则高斯-塞德尔迭代式(4.4)为

(4.6)

(4.7)

称为高斯-塞德尔迭代矩阵。

高斯-塞德尔迭代算法

高斯—塞德尔迭代的程序实现与雅可比迭代步骤大致相同,对于方程组 法中,假定雅可比迭代矩阵为 算法给出由 和 计算 表示 表示

,在前面的雅可比算 。下面的

,其迭代核心部分是
和对 初始化的部分。

的过程,省略了形成迭代矩阵

雅可比迭代的核心部分:

WHILE for(u:=1;u<=n,u++) x1[u]:=x2[u] for(i:=1;j<=n;j++) { s:=g[i]; for(j:=1;i<=n;i++)

{ s:=s+b[i]

x2[j]} //注意 x2[j]

}

ENDWHILE

在高斯-赛德尔迭代计算中并不需要形成迭代矩阵 在程序中令向量

,由式(4.5)可知在计算中只要形成矩阵

。它的核心部分是计算迭代式 ,计算中只需及时将 高斯-塞德尔迭代的核心部分: 放到 的位置上。

WHILE for(u:=1;u<=n;u++) x1[u]:=x2[u] for(i:=1;j<=n;j++) { s:=g[i]; for(j:=1;j<=n;j++)

{ s:=s+b[i] x2[i]:=s } ENDWHILE

x2[ j ]}

//注意 x2[j]

上列算法是在假定迭代收敛的前提下,使用当型(WHILE)结构控制循环。更一般地,可将上列算法 中 WHILE 循环改为 FOR 循环,通过控制循环次数和观测计算误差终止循环。届将上列算法中 WHILE 语 句改为

WHILE 这时在程序中要增加循环变量的设定和运算。

循环次数

判断高斯塞德尔迭代收敛的方法与判断雅可比迭代收敛类似,一方面从高斯-塞德尔迭代矩阵 信息,当 出判断。 定理 4.5 若方程组系数矩阵 A 为列或行对角优时,则高斯塞德尔迭代收敛。 或 的某种范数

获取

时,迭代收敛;另一方面,直接根据方程组系数矩阵的特点作

定理 4.6 若方程组系数矩阵 A 为对称正定阵,则高斯塞德尔迭代收敛。

例 4.3 方程组

中,

,

证明当

时 Gauss-Seidel 法收敛,而 Jacobi 迭代法只在

时才收敛。

解:对

法,因为

是对称矩阵,因此只要证



正定即可,由顺序主子



而 故 对称正定, 法收敛。



于是得到





法,可根据定理 4.2,由于迭代矩阵





法收敛的充要条件,故

法只在

时才收敛。



时,

法收敛,而



法不收敛,此时

不是正定的。

4.3
逐次超松弛迭代计算

逐次超松弛(SOR)迭代法

方程组 角矩阵。

的雅可比迭代形式



其中

是下三角矩阵,

是上三

得高斯-塞德尔的迭代形式:



,有

这样 加上修正量

可视为

的修正量,而

恰是由

加修正量

而得,如果将

改为

)

乘一个因子

,迭代格改为:

整理得

(4.8) 这里 为修正因子,称为松弛因子,而式(4.8)称为松弛迭代。迭代的分量形式为

(4.9) 式(4.9)称为松弛迭代的计算格式。 例 4.4 给定方程组

用 SOR 法求解,取 解:用 SOR 迭代公式可得

取初始值:

如果用高斯-赛德尔迭代法

迭代 72 次得:

用 SOR 迭代法

,只须迭代 25 次即得:

逐次超松弛迭代算法
下列算法假定迭代矩阵收敛,否则要在 WHILE 循环中增加判断条件。

1.定义和输入系数矩阵

与常数项向量

的元素,输入松弛因子

的值。

2.FOR i:=1,2,…,n

// 假定

,形成常数项向量

FOR

}

3.

4. WHILE

;

} ENDWHILE

5.输出

.

松弛迭代矩阵

将式(4.9)中含有



的项分别放在方程两边:



代入上式,有

则松弛因子为

的迭代矩阵为

其中

定理 4.7 逐次超松弛迭代收敛的必要条件



定理 4.8 若

为正定矩阵,则当

时,逐次超松弛迭代恒收敛。 究竟取多少为最佳, 这

以上定理给出了逐次超松弛迭代因子的范围。对于每个给定的方程组,确定 是比较困难的问题,对某些特定的方程组,我们可以得到一些理论结果。

通常,把 迭代称为松弛迭代。

的迭代称为亚松弛迭代,把

的迭代称为高斯-塞德尔迭代,而把



4.4

逆矩阵计算

在线性代数中逆矩阵是按其伴随矩阵定义的,若 的伴随矩阵。要计算 采用。通常对 个

则方阵

可逆,且

,其中



阶的列式才能得到一个伴随矩阵,在数值计算中因其计算工作量大而不被 化成 的过程中得到 。在数值计算中,这仍然是

做行的初等的效换,在将

一种行之有效的方法。

由逆矩阵的定义



,有

化为 个方程组

j

是第 个分量为 1,其余分量为 0 的 维向量。或记为:



用直接法或迭代法算出

也就完成了逆矩阵

计算。

如果依次对

用高斯若尔当消元法,组合起来看有(当然也能组合起来做):

这正是在线性代数中用初等变换计算逆矩阵的方法。 由此可见,计算一个 阶逆矩阵的工作量相当于解 个线性方程组。在数值计算中常常将计算矩阵逆 的问题转化为解线性方程组的问题。

例如, 已知方阵 的乘积得到 项解出 。

和向量 ; 而将

有迭代关系式

, 在计算中不是先算出

, 再作



作为线性方程组系数矩阵, 求解方程组

作为常驻数

4.5
程序 4.1 用雅可比迭代求解方程组:

程序示例

算法描述

输入系数矩阵

及常数项向量 c。

按雅可比迭代公式:

求解



程序源码
///////////////////////////////////////// // Purpose:雅可比迭代求解线性方程组 //

///////////////////////////////////////// #include<stdio.h> #include<stdio.h> #define MAX-N 20 #define MAXREPT 100 #define epsilon 0.00001 int main( ) { int n; int I, j, k; double err; static double a [MAX-N] [MAX-N],b[MAX-N] [MAX-N],c[MAX-N],g [MAX-N]; static double x[MAX-N],nx[MAX-N]; printf(〃\nInput n value(dim of Ax=c):〃); //输入方程的维数 scanf(〃%d〃,&n); if(n>MAX-N) //求解精度 //方程的最大维数

{ printf(〃The input n is larger then MAX-N,please redefine the MAX-N.\n〃); } if(n<=0) { printf(〃please input a number between 1 and %d.\n)〃, MAX-N};return 1;} //输入 Ax=c 的系数矩阵 A. PRINTF(〃Now input the matrix a(I,j,)=0,…,%d:\n〃,n-1);

for(i=0; i<n; i++) for(j=0; j<n; j++) {

//形成

迭代矩阵 b

b[i][j]=-a[i][j]/a[i][i]; g[i]=c[i]/a[i][j]; //为了简化程序,假设 a[i][i]!≠0,

//否则要附加对 a[i][i]的判断及其相应的处理 { for(i = 0; i< MAXREPT; i++) { for(j=0; j<n; j++ ) nx[j]=g[j]; for(j=0, j<n; j++ ) { for(k=0, k<n; k++ ) { if (j==k)continue; nx[j]+ =b[j][k]*x[k]; } //迭代

} err = 0 ; for(j=0, j<n; j++ ) if (err<fabs(nx[j]-x[j]))err=fabs(nx[j]-x[j]); for(j=0, j<n; j++ ) x[j]=nx[j]; if(err<epsilon) { printf(“Solve…x_I=\n”); //输出 //误差

for(i=0, i<n; i++ )printf(“%f\n”,x[i]); return 0; { } printf(“After% d repeat, no result…\ n”, MAXREPT ); return 1; { /输出

计算实例
用雅可比方法解下列方程组:

程序输入输出
Inprt n value(dim of Ax=c):3 Now input the matrix a(j, j), i, j =0,…2: 64 -13 -1 2 –90 1 1 1 40 Now input the matrix c(i), i= 0,…, 2:

14 -5 20 Solve…x_i= 0.229547 0.066130 0.492608 程序 4.2 用逐次超松弛迭代( 作参数)求解方程组:

算法描述

输入矩阵

及列向量 c。 的超松弛迭代公式:

按因子为

求解



程序源码
///////////////////////////////////////// // Purpose:超松弛迭代求解线性方程组 // ///////////////////////////////////////// # include<stdio.h> # include<math.h> # define MAX_N 20 / /方程的最大维数

# define MAXREPT 100 # define epsilon 0.00001 int main( ) { int n; int i, j, k ; double err, w; static double a[MAX_N][MAX_N], b[MAX_N][MAX_N],c[MAX_N],c[MAX_N],g[MAX_N] ; static double a[MAX_N] , nx[MAX_N]; PRINTF(“\nlnput n value(dim of Ax=c ):”); //输入方程的维数 Scanf(“% d”, &n); If (n>MAX_N) { printf(“The input n-is larger then MAX_N, please redefine the MAX_N.\n”); return 1; } if (n<=0) { printf(“Please input a number between 1 and % d. \n”, MAX_N); return 1;} / /求解精度

/ /输入

的 A 矩阵

printf(“Now input the matrix a(i, j), i, j =0,…% d: \n”, n-1); for(i=0; i<n; i++) for(j=0; j<n; j++) scanf(“%1f” , &s[i][j]); printf(“Now input the matrix c(i), i =0,…% d: \n”, n-1); for(i=0; i<n; i++) scanf (“%1f” , &c[i]);

printf(“Now input the w value:”); scanf(“%1f” , &w); if (w<1|| w>=2) } printf(“w must between 1and 2. \ n”); return 1; { for(i=0; i<n; i++) for(j=0; j<n; j++) { b[i][j]=-a[i][j]/a[i][i]; g[i]=c[i]/a[i][i]; { for(i=0; i<MAXREPT; i++) { for(j=0; j<n; j++) nx[j]=g[j]; for(j=0; j<n; j++) { for(k=0; k<j; k++) //改造 x {k+1}=bx_{k}+g 迭代矩阵

{ err=0;

for(j=0; j<n ; j++) if(err<fabs(nx[j]-x[j]))err=fabs(nx[j]-x[j]); for(j=0; j<n ; j++) x[j]=nx[j]; if(err<epsilon) { printf(“Solve … x_i = \ n”) ; //输出 //误差

for(i=0; i<n ; i++)printf(“%f\n”,x[i]); return 0; } { printf(“After % d repeat, no result … \ n”, MAXREPT); return 1 ; //输出

计算实例
解下列方程组:

程序输入输出
Input n value(dim of Ax=c): 3 Now input the matrix a(i, j), i, j=0,…,2: 64 -3 -1 2 -90 1 1 1 40 Now input the matrix c(i), i=0,…,2: 14 -5 20 Now input the w value: 1

Solve… x_i= 0.229547 0.066130 0.492608

第5章
5.1


引言



在实际问题中,有时只能给出函数 能给出 来逼近函数 的具体解析表达式, 或者

在平面上的一些离散点的值

,而不

的表达式过于复杂而难于运算。 这时我们需要用近似函数

,在数学上常用的函数逼近的方法有: 插值。 一致逼近。 均方逼近或称最小乘法。

? ? ?

什么是插值?简单地说,用给定的未知函数 与 在给定点的函数值相等,则称函数

的若干函数值的点构造

的近似函数

要求

为插值函数。例如:在服装店订做风衣时,选择好

风衣的样式后,服装师量出并记下你的胸围、衣长和袖长等几个尺寸,这几个尺寸就是风衣函数的插值点 数值,在衣料上画出的裁剪线就是服装师构造的插值函数 造合乎身材的插值函数。 ,裁剪水平的差别就在于量准插值点和构

定义 5.1

为定义在区间

上的函数, ,满足

,为



个互不相同的点, 为

给定的某一函数类。若 上有函数

则称



关于节点

在 上的插值函数;称点

为插值节点;称 称为被插函数。

,为插值型值点,简称型值点或插点;

这样,对函数

在区间

上的各种计算,就用对插值函数

的计算取而代之。

构造插值函数需要关心下列问题:
?

插值函数是否存在?

? ?

插值函数是否唯一? 如何表示插值函数?

如何估计被插函数与插值函数的误差?

5.2

拉格朗日(Lagrange)插值

可对插值函数

选择多种不同的函数类型,由于代数多项式具有简单和一些良好的特性,例如,

多项式是无穷光滑的,容易计算它的导数和积分,故常选用代数多项式作为插值函数。

5.2.1 线性插值

问题 5.1 给定两个插值点

其中

,怎样做通过这两点的一次插值函数?

过两点作一条直线,这条直线就是通过这两点的一次多项式插值函数,简称线性插值。如图 5.1 所示。

图 5.1 线性插值函数 在初等数学中,可用两点式、点斜式或截距式构造通过两点的一条直线。 下面先用待定系数法构造插值直线。

设直线方程为

,将

分别代入直线方程

得:



时,因

,所以方程组有解,而且解是唯一的。这也表明,平面上两个点,有且仅

有一条直线通过。用待定系数法构造插值多项式的方法简单直观,容易看到解的存在性和惟一性,但要解 一个方程组才能得到插值函数的系数,因工作量较大和不便向高阶推广,故这种构造方法通常不宜采用。



时,若用两点式表示这条直线,则有:

(5.1) 这种形式称为拉格朗日插值多项式。





称为插值基函数,计算



的值,易见

(5.2)

在拉格朗日插值多项式中可将 点的作用和地位都是平等的。

看做两条直线



的叠加,并可看到两个插值

拉格朗日插值多项式型式免除了解方程组的计算,易于向高次插值多项式型式推广。

线性插值误差

定理 5.1 记

为以 ,设

为插值点的插值函数, 一阶连续可导, 在

。这里 上存在,则对任意给定的 ,

至少存在一点

,使

(5.3)

证明 令

,因



的根,所以可设

对任何一个固定的点 ,引进辅助函数







由定义可得

,这样

至少有 3 个零点,不失一般性,假定 在每个区间至少存在一个零点,不妨记为 和

,分别在 ,即 ,

和 和 。

上应用洛尔定理,可知 ,对 在

上应用洛尔定理,得到



上至少有一个零点

现在对

求二次导数,其中

的线性函数),故有

代入 ,得

所以



5.2.2 二次插值

问题 5.2 给定三个插值点 (抛物线)插值多项式?



,其中

互不相等,怎样构造函数

的二次的

平面上的三个点能确定一条次曲线,如图 5.2 所示。

图 5.2 三个插值点的二次插值 仿造线性插值的拉格朗日插值,即用插值基函数的方法构造插值多项式。设

每个基函数

是一个二次函数,对

来说,要求

是它的零点,因此可设

同理



也相对应的形式,得



代入

,得

同理将

代入

得到



的值,以及



的表达式。

也容易验证:

插值基函数仍然满足:

二次插值函数误差:

上式证明完全类似于线性插值误差的证明,故省略。

插值作为函数逼近方法,常用来作函数的近似计算。当计算点落在插值点区间之内时叫做内插,否则 叫做外插。内插的效果一般优于外插。

例 5.1 给定 和 解:构造线性插值函数:

。 构造线性插值函数并用插值函数计算

分别将

代入上式,得

,准确值

,准确值

例 5.2 给定 。 解:

。构造二次插值函数并计算

,准确值

例 5.3 要制做三角函数的函数

值表,已知表值有四位小数,要求用线性插值引起的截断误差不

超过表值的舍入误差,试决定其最大允许步长。 解:设最大允许步长

5.2.3 n 次拉格朗日插值多项式

问题 5.3 给定平面上两个互不相同的插值点 给定平面上三个互不相同的插值点 面上 个互不相同的插值点

,有且仅有一条通过这两点的直线; ,有且仅有一条通过这三个点的二次曲线;给定平 ,互不相同是指 互不相等,是否有且仅有一条

不高于 次的插值多项式曲线,如果曲线存在,那么如何简单地作出这条 次插值多项式曲线?

分析: 次多项式 通过给定平面上 个互不相同的插值点

,它完全由

个系数 ,则

决定。若曲线 应满足

,事实上一个插值点就是一个插值条件。



依次代入

中得到线性方程组:

(5.4) 方程组的系数行列式是范德蒙(Vandermonde)行列式:

当 惟一。

互异时,

,所以方程组(5.4)的解存在且惟一。即问题 5.3 的解存在而且

通过求解(5.4)得到插值多项式

,因其计算量太大而不可取,仿照线性以及二次插值多项式

的拉格朗日形式,我们可构造 次拉格朗日插值多项式。

对于

个互不相同的插值节点 。

,由 次插值多项式的惟一性,可对每个插值节点

作出相应的 次插值基函数

要求



,的零点,因此可设





代入

,得到

(5.5)

作其组合:

(5.6)

那么

不高于 次且满足

,故

就是关于插值点 的拉格朗日基函数。

的插值多项式,这种插值形式称为拉格朗日插值多项式。

称为关于节点

例 5.4 给出下列插值节点数据,做三次拉格朗日插值多项式,并计算

(0.6)。

-2.00 17.00
解:拉格朗日插值基函数为:

0.00 1.00

1.00 2.00

2.00 17.00

三次拉格朗日插值多项式:

n 次插值多项式的误差

定理 5.2 设 当



上过

的 次插值多项式,

互不相等,

时,则插值多项式的误差:

其中

(5.7)

证明*:记 于是可设

。由于

,因而



的根,

下面的目标是算出

,为此引入变量为 的函数



(5.8)



,得



,由定义 , 由洛尔定理, 在



至少有

个零点, 由于

相邻的两个零点之间至少有一个零点, 即

至少有

个零点。同理再对 至少有一个零点 。

应用洛尔定理,即

至少有 个零点,反复应用洛尔定理得到

另一方面,对



阶导数,有



,有

得到

(5.9)

由于

的零点



的零点

有关,因而

为 的函数。

若|

可表示为

(5.10)

由(5.9)式可以看到,当

是不高于 次的多项式时,

,即



对于函数 拉格朗日基函数 满足

,关于节点

的拉格朗日插值多项式就是其本身,故



,得到



定理 5.2 给出了当被插函数充分光滑时的插值误差或称插值余项表达式,但是,在实际计算中,并不 知道 的具体表示,难以得到 的形式或较精确的界限 ,因此也难以得到界 。在实

际计算中,可对误差运用下面的事后估计方法。

给出

个插值节点 。在

,任选其中的

个插值节点,不妨取 个插值点,不妨取

,构造一 ,

个 次插值多项式,记为

个插值节点中另选 。由定理 2 可得到

构造一个 次插值多项式,记为

(5.11)

(5.12)



在插值区间内连续而且变化不大,有

,则

从而可得到

(5.13)

(5.14)

拉格朗日插值多项式的算法
下面用伪码描述拉格朗日插值多项式的算法。

1: 输入: 插值节点控制数 , 插值点序列

, 要计算的函数点 , 及变量



2:FOR i :=0,1,…,n {tmp:=1;

//i 控制拉格朗日基函数序列

2.1 FOR j:=0,1,…,i-1,i+1,…,n

{ //对于给定 ,计算拉格朗日基函数



}

// tmp 表示拉格朗日基函数

2.2 }

3:输出

的计算结果



5.3

埃特金逐步插值法

我们已经知道,拉格朗日插值多项式的余项为

其中

但在实际插值过程中,由于

一般无法计算,因而实际的误差也就无法得知。那么,如何根

据所要求的精度来选取插值点的多少呢?这就是埃特金(Aitken)逐步插值所要解决的问题。 首先介绍一个误差的事后估计法。



为通过

点的 次插值多项式;

为通过 ,而后者还通过

点的 ,即

次插值多项式。这两个 次插值多项式都通过

,其中前者还通过

它们所通过的插值线结点有 个是共同的,只有一个是不同的。若插值多项式用 一个下标 表示插值多项式通过 个插值结点

表示,则其中第 。

;第二个下标 表示还通过的插值结点

现在考虑插值多项式



,由于它们所通过的插值结点中有一个不同的

, 因

此其余项也不同。设它们的余项分别为

其中



虽然不等,但由于在插值区间内这两个插值多项式所通过的插值结点大 阶导数相差不多(假设 与 离得很近),可以近似看作相

部分相同,只有一个不同,所以这两个 等,即

由此可以得到这两个插值多项式的余项比为

从上式中解出





由上式可以看出, 次插值多项式

的误差可以由具有 个相同的结点与一个不同结点的两个 的近似

次插值多项式的差来估计。并且还可以看出,如果将这个误差加上,则可以得到更精确的 值。可以证明,这更精确的插值多项式为 次插值多项式,即

由上式可以看出,两个 次插值多项式 项式



经过线性组合后,便得到

次插值多

。这种方法称为埃特金逐步线性插值。它的优点在于可根据精度的要求逐步提高插值的阶,

在插值过程中,由于都是线性运算,因而计算也方便。其计算格式如图 5.3 所示。

图 5.3 埃特金逐步线性插值的计算格式

例 5.8 设函数 已知



的近似值。

解:根据埃特金逐步线性插值的计算格式,其计算结果如表 5.2 所示。 表 5.2

K xk

f(xk)

P0,k(x)

P1,k(x)

P2,k(x)

P3,k(x)

0 0.3 0.29854 1 0.4 0.39646 0.457195 2 0.5 0.49311 0.456134 0.456537 3 0.6 0.58813 0.454900 0.456484 4 0.7 0.68122 0.453502 0.456432 0.456557 0.456557 0.456557

由表 5.2 可知,若取

,则

所以有

此时,

实际上已用不着计算了。

若取

,则

所以有 此时,以后的值也不必计算了。

埃特金逐步线性插值的算法

输入:插值结点

;插值点 及精度要求 。

输出:插值点 处的函数值



在上述算法中,

分别表示

。 当某相邻两个



差的绝对值小于 时,插值过程就停止。如果所有的插值结点都已用完,但还未满足精度要求,插值过程 也停止。

5.4

埃尔米特(Hermite)插值

在构造插值时,如果不仅要求插值多项式节点的函数值与被插函数的函数值相同,还要求在节点处的 插值函数与被插函数的一阶导数的值也相同,这样的插值称为埃尔米特插值或称密切插值。

常用埃尔米特插值描述如下:设 次的多项式函数 满足

具有一阶连续导数,以及插值点

,若有至多为

则称



关于节点

的埃尔米特插值多项式。

问题 5.7 5.7:给定 数值和一阶导数值的埃尔米特插值多项式?

;怎样构造给定两个节点的函

分析:用 4 个条件,至多可确定三次多项式。设满足插值条件的三次埃尔米特插值多项式为:

将插值条件代入

得到线性方程组:

因为方程组的系数行列式

所以方程组有解,可惟一解出

,即关于节点

的埃尔米特插值多项式存在惟一。类似 。

于拉格朗日插值多项式样构造手法,也可通过插值基函数作出





可设

同理要











由上述要求,对

来说,至多为三次多项式,

是它的二重根,可设

(5.20)

利用

解出

同理可得



,可设





,算出

。所以

(5.21) 同理

综上所述,得到以

为节点的埃尔米特插值为

(5.22)

容易证明,当

时插值误差为

(5.23)

如果要构造 方法类似。

关于

个节点



米特插值多项式,手法与构造两个节点的

这里

分别为不高于

次插值多项式,分别满足

及 由此可得到

(5.24)

这里

为关于节点

的拉格朗日基函数。

容易证明,当

时,误差为:

(5.25)

例 5.8 给定

,求埃尔米特插值多项式,并计算



解:

显然本题不必计算



利用构造基函数方法做插值多项式被广泛地应用在不同的插值条件中。

例 5.9 给定

构造二次插值多项式函数。

解:设

这里

均为不高于二次的多项式,它们分别满足

于是

可表示为



,得



,得

所以



所以

同理



所以

所以

用牛顿差商插值也能构造埃尔米特插值。对给定的插值点的函数值和一阶导数值 定义序列 即

计算一阶差商时:





,取

即构造差商表中用

代替

其余差商计算公式不变,得到差商型埃尔米特插值公式:

(5.26)

其中





例 5.10 用下列数据构造埃尔米特插值多项式,并计算 f(1.36)。 1.2 0.6 0.5 1.4 0.9 0.7 1.6 1.1 0.6

解:计算差商。 1.2 1.2 0.6 0.6 0.5

1.4 1.4 1.6 1.6

0.9 0.9 1.1 1.1

1.5 0.7 1.0 0.6

5 -4 1.5 -2 -45 13 -17.5 145 -76.25 -553.125

例 5.11 求

使



插值多项式及其余项表达式。

解:这里给出了四个条件故可构造三次插值多项式 Newton 均差插值,令

,可用

显然它满足条件

, 为待定参数。



可得

解得: 于是就可得到插值多项式。 它的余项为





之间,而



5.5

分段插值

5.5.1 龙格(Runge)现象
在构造插值多项式时,根据误差表达式(5.9),你是否认为多取插值点总比少取插值点的效果好呢? 答案是不一定。如果被插函数是高次多项式,那么多取插值点比少取插值点效果好;但对有些函数来说, 有时点取的越多,效果越不近人意。请看下面的例子。

给定函数



。构造 10 次插值多项式



对 多项式

作等距分割,取

,构造 对

,10 次插值 的逼近效果较好, 在 =-0.90,

如图 5.4 所示。 从图中可以看到, 在零点附近,

-0.70,0.70,0.90 时误差较大。

图 5.4



下面列出



的几个插值点的数值: -0.90 -0.70 -0.50 -0.30

1.57872 0.4706

0.07547

0.13793

0.30769 0.23535

-0.22620 0.25376

这个例子是由龙格(C.Runge)提出的,也称插值多项式在插值区间发生剧烈振荡的现象为龙格现象。 龙格现象揭示了插值多项式的缺陷。它说明高次多项式的插值效果并不一定优于低次多项式的插值效果。

在插值过程中,误差由截断误差和舍入误差组成。式(5.9)给出的是截断误差,它是插值函数 原函数 的误差。另外由节点



计算产生的舍入误差,在插值计算过程中可能被扩散或放大,这就是

插值的稳定性问题。而高次多项式的稳定性一般比较差,这从另一角度说明了高次插值多项式的缺陷。

5.5.2

分段线性插值

既然增加插值点并不能提高插值函数的逼近效果,那么采用分段插值的效果又如何呢?

若对给定区间

作分割 ,则

, 在每个小区间

上作



为节

点的线性插值,记这个插值函数为

(5.27)

把每个小区间的线性插值函数连接起来,就得到了 分段线性函数 。 在

以剖分(节点) 是平面上以点



上为一个不高于一次的多项式,事实上 时有

为折点的折线。由线性插值误差公式,当

(5.28)

因而

其中



于是,当区间分割加密, 续,分段线性插值序列就能收敛于 。

时,分段线性插值收敛于

。事实上,只要



分段线性插值算法简单,只要区间充分小,就能保证它的误差要求。它的一个显著优点是它的局部性 质,如果修改了某节点 的缺点是在插值节点处不光滑。 的值,仅在相邻的两个区间 受到影响。分段线性插值

图 5.5 给出分段线性插值 于整体的拉格朗日插值的效果。

(虚线表示)和

的图形,可以看到分段线性插值的效果明显好

图 5.5 分段线性插值



例 5.12 对下列数据作分段线性插值,并计算

(1.2 ), 3 6 9

(3.3 )。

-3 -1 2 12 12 1

12

解: ∵ 1.2 [-1,2]



(1.2 ) =

=

×5+

×1=2.0667

∵ 3.3 [3,9]



(3.3 ) =

=

×6+

× 12 = 6.3

5.6

三次样条函数

在制造船体和汽车外形等工艺中传统的设计方法是,首先由设计人员按外形要求,给出外形曲线的一 组离散点值 将压铁放在点 ,施工人员准备好有弹性的样条(一般用竹条或有弹性的钢条)和压铁, 的位置上,调整竹条的形状,使其自然光滑,这时竹条表示一条插值曲线,我们称

为样条函数。从数学上看,这一条近似于分段的三次多项式,在节点处具有一阶和二阶连续微商。样条函 数的主要优点是它的光滑程度较高,保证了插值函数二阶导数的连续性,对于三阶导数的间断,人类的眼 睛已难以辨认了。样条函数是一种隐式格式,最后需要解一个方程组,它的工作量大于多项式拉格朗日型 式或牛顿型式等显式插值方法。

定义 给定区间

上 。若

个节点 满足 上有连续的二阶导数,则称 为样条节点。

和这些点上的函数值, ; 为 在每个小区间 关于剖分 上至多是一

个三次多项式;



的三次样条插值函数,称

要在每个子区间

上构造三次多项式

, 共需要 条件;用每个内点的关系建立条件

个条件, 由插值条件

, 提供了



又得到

个条件;再附加两个边界条件,即可惟一确定样条函数了。用待定系数法确定了构造样

条函数的存在性和惟一性。在具体构造样条函数时一般都不使用计算量大的待定系数法。下面给出构造三 次样插值的 关系式和 关系式的方法。

5.6.1

三次样条插值的 M 关系式

引入记号 系式,用一阶导数表示样条插值函数时称为小

。 用节点处二阶导数表示样条插值函数时称为大 关系式。



问题 5.8 给定插值点 关系式?

, 怎样构造用二阶导数表示的样条插值函数, 即怎样构造

假设 线性插值函数:

。由于



上为线性函数,故在

上做

的分段



,得到

(5.29)



积分两次有

(5.30)



代入式(5.27)可解出

故在

上有

(5.31)

在每个小区间上具有不同的表达式,但由于

在整个区间

上是二阶光滑的,故有

列出每一个关系式

,再经计算得:

(5.32) 其中:

由式(5.32)得到

个未知数的

个方程组。现补充两个边界条件,使方程组只有惟一解。下面

分三种情况讨论边界条件。

(1)给定 知量

的值 ,即

时,称为自然边界条件),此时

阶方程组有

个未

(5.33)

(2)给定 到另外两个方程:

的值,它们分别代入



中的表达式,得

于是需要解

阶的方程组:

(5.34)

(3) 被插函数以 即

为基本周期时, 即

, 即



。此时化为 个变量、 个方程的方程组。

样条插值构造的

关系式是对角占优的三对角带状矩阵,可用第 3 章中的追赶法求解。

例 5.13 给出离散数值表: 1.1 1.2 1.4 1.5 1.8000

0.4000 0.8000 1.6500



,构造三次样条插值的

关系式,并计算 f(1.25)。

解:由题中

的数值,计算得



的边界条件,得

解得



因此,三次样条插值的分段表达式为

特别地, 详细的程序和算例请看 5.8 节。



5.6.2

三次样条插值的 m 关系式

问题 5.9 给定插值点 构造 关系式?

,怎样构造用节点处一阶导数表示的样条插值函数,即怎样

对给定的插值点 特插值,那么在整个

,先假定已知 上是分段的埃尔米特插值,在 上

,在每个小区间 的表达式为

上做埃尔米

通过

得到方程组

其中:

再附加两个边界条件,即可解出

的值。附加的边界条件情况同

关系式中的讨论类似,不再详述。



作样条插值,插值效果见图 5.6。可以看到样条插值效果优于分段插值效果。

图 5.6 样条插值图示

5.7

程序示例

程序 5.1 给定

,构造牛顿插值多项式

互不相同。

算法描述

输入 值,及(

;记



计算差商

其中 对给定的 ,由

计算

的值;

输出



程序源码
////////////////////////////////////////////// // purpose:(x_i,y_i)的牛顿插值多项式 // ///////////////////////////////////////////// # include <stdio.h>

# define MAX_N 20 typedef struct tagPOINT { double x; double y; } POINT; int main ( ) { int n; int i,j;

// 定义(x_i,y_i)的最大维数 // 点的结构

POINT points [MAX_N + 1];double diff [MAX_N + 1]; double x,tmp,newton = 0; printf (“\ nInput n value:”); // 输入被插值点的数目 scanf (“% d”,&n); if (n>MAX_N ) { printf (“The input n is larger then MAX_N,please redefine the MAX_N.\ n”); return 1; } if (n<= 0) { printf (“Please input a number between 1 and % d. \ n”,MAX_N ); return 1; } // 输入被插值点(x_i,y_i) printf (“Now input the(x_i,y_i),i=0,…,% d: \ n”,n); for (i=0;i<= n;i++)

scanf (“% lf % lf”,&points[i].x, &points[i].y); printf (“Now input the x value:”); // 输入计算牛顿插值多项式的 x 值 scanf (“%lf”,& x); for (i=0;i<= n;i++) diff [i] = points [i].y; for (i=1;i<= n;i++) { for ( j = n;j>=i;j--) { diff [j] = ( diff [j]-diff [j-1]) / (points [j]. x – points [j-1-i].x); } } tmp = 1;newton = diff [0]; for (i=0;i<= n;i + +) { tmp = tmp * (x-points[i].x); newton = newton + tmp * diff [i+1]; } printf (“newton (% f ) = %f \ n”,x,newton); // 输出 return 0; } // 计算 f (x_0,…,x_n)的差商

计算实例
给定 sin11°= 0.190809,sin12°= 0.207912,sin13°= 0.224951,构造牛顿插值函数并计算 sin11°30′。

程序输入输出
input n value:2

Now input the(x_i,y_i),i=0,…,2: 11 0.190809 12 0.207912 13 0.224951 Now Input the x value:11.5 newton (11.500000) = 0.199369

程序 5.2 给定插值点 插值多项式 ,求在给定点 处

和二阶导数的端点值 的值。

,用

关系式构造三次样条

算法描述

1.输入 值,及 2.for i =1 to n-1



,要计算的函数点 ;

计算







其中 3.求解方程组



对给定的 ,由

4.计算出

的值;

5.输出



程序源码
///////////////////////////////////////////////////

// purpose:给定



值的三次样条插值多项式 //

/////////////////////////////////////////////////// # include <stdio.h> # define MAX_N 20 typedef struct tagPOINT { double x;double y; } POINT; int main ( ) { int n; int i,k; POINT points[MAX_N +1]; double h[MAX_N +1],b[MAX_N +1],c[MAX_N +1],d[MAX_N +1],M[MAX_N +1]; double u[MAX_N +1],v[MAX_N +1],y[MAX_N +1]; double x,p,q,S; printf (“\ nInput n value:”); // 输入插值点的数目 // 定义(x_i,y_i)的最大的维数 // 点的结构

scanf (“% d”,&n); if (n>MAX_N) { printf (“The input n is larger than MAX_N, please redefine the MAX_N. \ n”); return 1; } if (n<=0) { printf (“Please input a number between 1 and % d. \ n”,MAX_N); return 1 } // 输入插值点(x_i,y_i),M0 值和 Mn 值 printf (“Now input the (x_i,y_i),i=0,…,% d:\ n”,n); for (i=0;i<n;i++) scanf (“% lf % lf ”,&points[i].y); printf (“Now input the M[0] value:”); scanf (“% lf ”,&M[0]); printf (“Now input the M[n] value:”); scanf (“%lf ”,&M[n]); printf (“Now input the x value:”); //输入计算三次样条插值函数的 x 值 scanf (“%lf”,&x); if (x>points[n]. x | | x<points[0].x= { printf (“Please input a number between % f and % f. \ n”,points[0].x,points[0].x); return 1; } // 计算 M 关系式中各参数的值

h[0]=points[1].x-points[0].x; for (i=1;i<;i++ ) { h[i]=points[i+1].x-points[i].x; b[i]=h[i]/(h[i]+h[i-1]; c[i]=1-b[i]; d[i]=6*((points[i+1].y-points[i].y)/h[i] -(points[i].y-points[i-1].y)/h[i-1])/(h[i]+h[i-1]); } // 用追赶法计算 Mi,i=1,…,n-1 d[1]-= c[1]*M[0]; d[n-1]-=b[n-1]*M[n]; b[n-1]=0;c[1]=0;v[0]=0; for ( i=1;i<n;i++) { u[i]=2-c[i]*v[i-1]; v[i]=b[i]/u[i]; y[i]= (d[i]-c[i]*y[i-1])/u[i]; } for (i=1;i<n;i++) { M[n-i]=y[n-i]-v[n-i]* M[n-i+1]; } // 计算三次样条插值函数在 x 处的值 k=0; while (x>= points[k].x) k++; k=k-1

p=points[k+1].x-x;q=x-points[k].x; S=(p* p* p* M[k] +q* q* q* M[k+1]) / (6* h[k]) +( p* points [k].y+q* points[k+1].y) / h[k]-h[k]*(p* M[k] +q* M[k]+q* M[k+1])/6; printf (“S(%f ) = % f\ n”,x,S); // 输出 getchar ( ); return 0; }

计算实例

给定离散点(1.1,0.4),(1.2,0.8),(1.4,1.65),(1.5,1.8), 构造三次样条插值多项式 ,计算 。

,用

关系式

程序输入输出
Input n value:3 Now input the (x_i,y_i),i=0,…,3: 1.1 0.4 1.2 0.8 1.4 1.65 1.5 1.8 Now Input the M0 value:0 Now Input the Mn value:0 Now Input the x value:1.25 S (1.25000) = 1.033594

第6章

曲线拟合的最小二乘法
6.1 拟合曲线

通过观察或测量得到一组离散数据序列 数 逼近客观存在的函数

, 当所得数据比较准确时, 可构造插值函

,构造的原则是要求插值函数通过这些数据点,即 与 是相等的。

。此时,序列

如果数据序列

,含有不可避免的误差(或称“噪音”),如图 6.1 所示;如果数据 最优地靠近样点,即向 与 之间误差最小原则

序列无法同时满足某特定函数,如图 6.2 所示,那么,只能要求所做逼近函数 量 与 的误差或距离最小。按

作为“最优”标准构造的逼近函数,称为拟合函数。

图 6.1 含有“噪声”的数据

图 6.2 一条直线公路与多个景点 插值和拟合是构造逼近函数的两种方法。插值的目标是要插值函数尽量靠近离散点;拟合的目标是要 离散点尽量靠近拟合函数。

向量



之间的误差或距离有各种不同的定义方法。例如:

用各点误差绝对值的和表示:

用各点误差按模的最大值表示:

用各点误差的平方和表示:



(6.1)

其中

称为均方误差,由于计算均方误差的最小值的方法容易实现而被广泛采用。按均方误差达到极

小构造拟合曲线的方法称为最小二乘法。本章主要讲述用最小二乘法构造拟合曲线的方法。 在运筹学、统计学、逼近论和控制论中,最小二乘法都是很重要的求解方法。例如,它是统计学中估 计回归参数的最基本方法。 关于最小二乘法的发明权,在数学史的研究中尚未定论。有材料表明高斯和勒让德分别独立地提出这 种方法。勒让德是在 1805 年第一次公开发表关于最小二乘法的论文,这时高斯指出,他早在 1795 年之前 就使用了这种方法。但数学史研究者只找到了高斯约在 1803 年之前使用了这种方法的证据。 在实际问题中, 怎样由测量的数据设计和确定“最贴近”的拟合曲线?关键在选择适当的拟合曲线类型, 有时根据专业知识和工作经验即可确定拟合曲线类型;在对拟合曲线一无所知的情况下,不妨先绘制数据 的粗略图形,或许从中观测出拟合曲线的类型;更一般地,对数据进行多种曲线类型的拟合,并计算均方 误差,用数学实验的方法找出在最小二乘法意义下的误差最小的拟合函数。 例如,某风景区要在已有的景点之间修一条规格较高的主干路,景点与主干路之间由各具特色的支路 联接。设景点的坐标为点列 条直线。通过计算均方误差 ;设主干路为一条直线 最小值而确定直线方程(见图 6.2)。 ,即拟合函数是一

6.2
线性拟合

线性拟合和二次拟合函数

给定一组数据

,做拟合直线

,均方误差为

(6.2)

是二元函数,

的极小值要满足

整理得到拟合曲线满足的方程:

(6.3)

或 称式(6.3)为拟合曲线的法方程。用消元法或克莱姆法则解出方程:

a=

=

例 6.1 下表为 P. Sale 及 R. Dybdall 在某处作的鱼类抽样调查,表中 为鱼的数量, 用线性函数拟合鱼的数量和种类的函数关系。 13 11 40 13 15 10 42 14 16 11 55 22 21 12 60 14 22 12 62 21 23 13 64 21 25 13 70 24 29 12 72 17 30 14 23 31 16 34 36 17

为鱼的种类。 请

100 130

解:设拟合直线

,并计算得下表: 编 号 1 2 3 4 5

x
13 15 16 21 22

y
11 10 11 12 12

xy
143 150 176 252 264

x2
169 225 256 441 484

21 ∑

130 956

34 344

4420

16900

18913 61640

将数据代入法方程组(6.3)中,得到:

解方程得: = 8.2084, = 0.1795

拟合直线为:

= 8.2084 + 0.1795

二次拟合函数

给定数据序列

,用二次多项式函数拟合这组数据。



,作出拟合函数与数据序列的均方误差:

(6.4)

由多元函数的极值原理,

的极小值满足

整理得二次多项式函数拟合的法方程:

(6.5)

解此方程得到在均方误差最小意义下的拟合函数 方程的系数矩阵是对称的。当拟保多项式阶 或一些特殊算法以保护解的准确性。

。方程组(6.5)称为多项式拟合的法方程,法

时,法方程的系数矩阵是病态的,在计算中要用双精度

例 6.2 给定一组数据,如下表。用二次多项式函数拟合的这组数据。 -3 4 -2 2 -1 3 0 0 1 -1 2 -2 3 -5

解:设

,由计算得下表:

-3 -2 -1 0 1 2 3

4 2 3 0 1 2 5 1

-12 -4 -3 0 -1 -4 -15 -39

9 4 1 0 1 4 9 28

36 8 3 0 -1 -8 -45 -7

-27 -8 -1 0 1 8 27 0

81 16 1 0 1 16 81 96

将数据代入式(6.5),相应的法方程为:

解方程得:

=0.66667,

=-1.39286,

=-0.13095



= 0.66667-1.39286 -0.13095

拟合曲线的均方误差: 结果见图 6.3。

=3.09524

6.3

解矛盾方程组

在 6.2 节中用最小二乘法构造拟合函数,本节中用最小二乘法求解线性矛盾方程组的方法构造拟合函 数。

给定数据序列 些点,那么有矛盾方程组:

,作拟合多项式

,如果要求

过这

即:

(6.6) 我们作一辅助函数

这是自变量为 一个自变量的偏导数为 0。

的多元函数,要使

达到最小值,采用多元函数求权值的方法,

对每

整理成以

为未知数的线性方程组

整理成矩阵形式

,其中:

这是一个 盾方程组的最小二乘解。

的对称方程组,称为法方程。只要

非奇异,就可以得出唯一解。这就是矛

有什么快捷的方法来求法方程的解?

把矛盾方程组(6.6)写成矩阵形式

,其中

容易验证

,法方程就是



例如,拟合直线

的矛盾方程组

的形式如下:

化简得到与(6.3)相同的法方程:

在线性代数中,我们知道,关于方程组

,若秩



,则方程组无解,这时方程组称 极小意义下的解, 也就是在 的解就是矛盾方程组 在

为矛盾方程组。 在数值代数中对矛盾方程计算的是在均方误差 最小二乘法意义下的矛盾方程的解。 定理 6.1 将证明, 方程组 最小二乘法意义下的解。 定理 6.1

1.



行 列的矩陈,

为列向量,秩



称为矛盾方程的



法方程,法方程恒有解。

2.



的解,当且仅当

满足

,即

是法方程的解。

证明:

1)对

作行初等变换

,使



∴ 秩

而秩





∴ 方程组

有解而且解惟一存在。

2)设

满足

,任取

,则

由于 数据

是任取的,故法方程组

的解为极小问题

的解。事实上,对离散

作 次多项式曲线拟合,要计算

的极小问题。这与解矛盾方程组





的极小问题是一回事。在这里

故对离散数据 得:

所作的 次拟合曲线

,可通过解下列方程组求

例 6.3 给定一组数据,见下表,用二次多项式函数拟合这组数据。 -3 -2 -1 4 2 3 0 0 1 2 3

-1 -2 -5

解:记二次拟合曲线为

,形成法方程:



=

=

=

= 得到:

=

解方程得:

= 0.66667,

= -1.39286,

= -0.13095



= 0.66667-1.39286

-0.13095

例 6.4 给出一组数据,见下表。用最小二乘法求形如

的经验公式。 -1 4.7 2 4

x y

-3 14.3

-2 8.3

8.3 22.7

解:列出法方程:

=

=

= 法方程为:

解方程得到: =10.675, = 0.137

拟合曲线为:

10.678+0.137

有些非线性函数经过转换以后可化为线性函数计算。例如, 曲线为: 。



,则化拟合

例 6.5 求一个形如

的经验函数公式,使它能够和下列数据相拟合。 1 15.3 2 3 4 5 6 7 8

20.5 27.4 36.6 49.1 65.6

87.8 117.6

解:化经验公式为线性形式,对经验公式的两边取自然对数有

由矛盾方程组

得到法方程组



解方程组得:

∴ 拟合曲线的均方误差为:

计算结果见图 6.4。

图 6.4 拟合曲线图示 例 6.6 解矛盾方程组

解:写出法方程组

,即



解法方程得:

1.5917,

0.5899,

0.7572。

图 6.3 拟合曲线与数据序

6.4



有的实际问题中,已知数据 不同,在矛盾方程组中,一组 把每对 重大些。

不一定都是一次观测的结果,对于不同的

可能观测次数

确定一个方程,而最小二乘解对每个方程来说都还存在误差,这样,

都同等看待就不太合理,希望观测次数多的(即可靠性大些的)数据在矛盾方程组中占的比

为了使问题的提法更有一般性,通常把最小二乘法中偏差平方和改为加权平方和,使

最小,其中

称作权。

同样,可由误差函数

求极值,得到法方程

,其中



事实上,只需对矛盾方程组

作一些修改,每个方程乘上

,得到

,其中

这样,就使 小二乘解。

,仍然可用前面的方法,对

,左乘



求得最

6.5

用正交多项式作最小二乘拟合

我们不仅可用多项式来拟合函数,还可用一般的函数来拟合。

定义 6.1 若

,如果

当且仅当 为 上线性无关族。

时成立,则称



上线性无关,称

最小二乘法的一般提法为:对给定的一组数据 中找一个函数 ,使加权平地均和

,要求在函数类

,其中

这里 表示该点数据的重要程度。

是线性无关的函数族,



上的权函数,点

处的

求误差函数 值的必要条件

的极小值点

, 由多元函数极

得:

这是

个未知数

个方程的方程组,称为法方程式。

定义 6.2 6.2:

称为



关于点集

的内积。

这样,法方程式可简写为

,记为

,其中

称为克莱姆行列式,记作



定理 6.2

的充要条件是

线性无关。

证明:若存在

使

对此式两边分别取与

的内积得:

这是一个以 是

为未知数的齐次方程组,有非零解的充要条件是系数矩阵行列式等于零,于 全为 0,所以 线性无关。证毕。

的充要条件是方程有全零解,即

由于法方程有惟一解的充要条件是 惟一解的充要条件。特别当 函数族,因而必有最小二乘解。

,因而 取为 时,由于

线性无关也是法方程 是在



中的线性无关



上的多项式拟合,需要解一个

的线性方程组,且当 取得大一此时,法方程

的系数矩阵会出现病态。从系数矩阵 B 的形式看,里面的元素都是些内积,是否能取某些函数族,使对非 对角元素全变为 0?如果有这样的函数族,那么方程容易解,病态也得到改善。

定义 6.3 函数族

如果在点集

上满足



为点集

带权

的正交函数族。

例 6.7 三角函数族

上是正交函数族(权

).

实际上

,而

如果拟合函数在

上取,且

是正交函数族,则法方程式成为:

直接可得到

,不用解线性方程组了。



容易估算,是否病态也容易判断。

完成以上工作的关键在于如何构造正交函数族。 正交多项式是最简单的正交函数族。常用的正交多项式如:Chebyshev(切比雪夫)多项式、Legendre (勒让德)多项式等。

现在我们根据给定结点

及权函数

,造出带权

正交的多项式族

,用递推公式表示如下:

其中

这样给出的

是正交的,这一点可以用归纳法证明。

例 6.8 已知函数表,利用正交多项式求拟合多项。

1 2 3 4 10 18 1 1 1

4 26 1

解:设

所以:

所以:

所以:

得:

6.6

程序示例

程序 6.1 用线性函数

拟合给定数据



算法描述

输入

值,及



解方程组:

输出



程序源码
//////////////////////////////////////////// // purpose:(x_i,y_i)线性拟合函数 //////////////////////////////////////////// # include <stdio.h> # define MAX_N 25 typedef struct tagPOINT { double x; double y; } POINT; int main ( ) { int m; int i; POINT points [MAX_N ]; static double u11,u12,u21,u22,c1,c2; double a,b,tmp; printf (“\ nInput n value:”); // 输入点数 m scanf (“% d”,&m); if (m>MAX_N ) { printf (“The input n is larger then MAX_N,please redefine the MAX_N.\ n”); return 1; } if (m<= 0) //(x_i,y_i)的最大维数 // 点的结构 //

{ printf (“Please input a number between 1 and % d. \ n”,MAX_N ); return 1;} printf (“Now input the(x_i,y_i),i=0,…,% d: \ n”,m-1); // 输入点 for (i=0;i<= m;i++) { scanf (“% lf”,&tmp); printf [i].x=tmp; scanf (“%lf ”,& tmp); printf [i].y=tmp; // scanf (“%lf %lf ”,& printf [i].x,& printf [i].y); } for (i=0;i<= n;i++) { u21+ = points [i].x; // 列出方程 U(a,b)T = c

u22+ = points [i].x points [i].x; c1+ = points [i].y;

c2+ = points [i].y points [i].y; } u12 = u21; u11 = m; // 求解

a = ( cl u22-c2 u12) / (u11

u22-u12 u21);

b = ( cl *u21-c2 u11) / (u21

u12-u22 u11); // 输出

printf(“Solve:p(x) = %f +%f x \ n”,a,b); return 0;

} // ----------End of File---------

计算实例
用线性函数拟合下列数据(例 6.1)。 13 11 40 13 15 10 42 14 16 11 55 22 21 12 60 14 22 12 62 21 23 13 64 21 25 13 70 24 29 12 30 14 31 16 36 17

72 100 130 17 23 34

程序输入输出
input m value:21 Now input the (x_i,y_i),i=0,…,20: 13 11 15 10 16 11 21 12 22 12 23 13 25 13 29 12 30 14 31 16 36 17 40 13 42 14 55 22 60 14 62 21 64 21 70 24 72 17 100 23 130 34 Solve:p(x) = 8.208408 + 0.179522 x

程序 6.2 用形如

的函数拟合定数据



算法描述

输入

值,及



解方程组:

输出所求拟合函数



程序源码
///////////////////////////////////////////////////////// // purpose:(x_i,y_i)形如 y=aExp(bx)的拟合函数 ///////////////////////////////////////////////////////// # include <stdio.h> # include <math.h> # define MAX_N 25 typedef struct tagPOINT { double x; double y; } POINT; int main ( ) { int m; int i; POINT points [MAX_N ]; static double u11,u12,u21,u22,c1,c2; double A,B,tmp; printf (“\ nInput n value:”); scanf (“% d”,&m); if (m>MAX_N ) { printf (“The input n is larger then MAX_N,please redefine the MAX_N.\ n”); // 输入点数 m //(x_i,y_i)的最大维数 // 点的结构 //

return 1; } if (m<= 0) { printf (“Please input a number between 1 and % d. \ n”,MAX_N ); return 1;} printf (“Now input the(x_i,y_i),i=0,…,% d: \ n”,m-1); // 输入点 for (i=0;i<= m;i++) { scanf (“% lf ”,&tmp); printf [i].x = tmp; scanf (“%lf ”,& tmp); printf [i].y = tmp; // } for (i=0;i<= n;i++) { // 列出方程 X(a,b)T = Y scanf (“%lf %lf ”,& printf [i].x,& printf [i].y);

// 其中 Y 为对原方程取对数后的数据

u21+ = points [i].x;

u22+ = points [i].x points [i].x; c1+ = log (points [i].)y;

c2+ = points [i].y log ( points [i].y); } u12 = u21; u11 = m;

// 求解

a = ( cl u22-c2 u12) / (u11

u22-u12 u21);

b = ( cl u21-c2 u11) / (u21

u12-u22 u11);

printf(“Solve:p(x) = % lfexp (%lf x) \ n”,exp(A),B); // 输出 return 0; } // ----------End of File---------

计算实例

求一个形如

的经验函数公式,使它能够和下列数据相拟合(例 6.5)。 1 2 3 4 5 6 7 8

15.3 20.5 27.4 36.6 49.1 65.6 87.8 117.6

程序输入输出
input m value:8 Now input the(x_i,y_i),i=0,…,7: 1 15.3 2 20.5 3 27.4 4 36.6 5 49.1 6 65.6 7 87.8 8 117.6 Solve:p(x) = 11.437069exp(0.291215x)

第7章

数值微分和数值积分
7.1 数值微分

7.1.1 差商与数值微分

当函数

是以离散点列给出时,当函数的表达式过于复杂时,常用数值微分近似计算

的导数

。在微积分中,导数表示函数在某点上的瞬时变化率,它是平均变化率的极限;在几何上可解释为 曲线的斜率;在物理上可解释为物体变化的速率。

以下是导数

的三种定义形式:

(7.1) 在微积分中,用差商的极限定义导数;在数值计算中返璞归真,导数取用差商(平均变化率)作为其 近似值。 最简单的计算数值微分的方法是用函数的差商近似函数的导数, 即取极限的近似值。 下面是与式 (7.1) 相应的三种差商形式的数值微分公式以及相应的截断误差。

向前差商
用向前差商(平均变化率)近似导数有:

(7.2)

其中

的位置在

的前面,因此称为向前差商。同理可得向后差商、中心差商的定义。

由泰勒展开

得向前差商的截断误差:

向后差商
用向后差商近似导数有:

(7.3) 与计算向前差商的方法类似,由泰勒展开得向后差商的截断误差:

中心差商
用中心差商(平均变化率)近似导数有:

(7.4) 由泰勒展开

得中心差商的截断误差:

差商的几何意义

微积分中的极限定义

,表示



处切线的斜率,即图 7.1 中

直线 一条过

的斜率;差商

表示过



两点直线

的斜率,是

的割线。可见数值微分是用近似值内接弦的斜率代替准确值切线的斜率。

图 7.1 微商与差商示意图

例 7.1 给出下列数据,计算 -0.02 0.04 0.06 5.06 0.8

, 0.10 5.055

5.07 5.065 5.05

解:

(5.07-5.06)/(0.04-0.02)= 0.5

(5.05-5.07)/(0.08-0.04)= -0.5

(5.05-5.055)/(0.08-0.10)= 0.25



(0.10) -

(0.06))/(0.10-0.06)= 18.75

设定最佳步长
在计算数值导数时,它的误差由截断误差和舍入差两部分组成。用差商或插值公式近似导数产生截断 误差,由原始值 的数值近似产生舍入误差。在差商计算中,从截断误差的逼近值的角度看, 越小,

则误差也越小;但是太小的 最小呢?

会带来较大的舍入误差。怎样选择最佳步长,使截断误差与舍入误差之和

一般对计算导数的近似公式进行分析可得到误差的表示式,以中心差商为例,截断误差不超过

而舍入误差可用量

估计(证明略),其中 是函数

的原始值的绝对误差限,总误差为



时,总误差达到最小值,即

(*) 可以看到用误差的表达式确定步长,难度较大,难以实际操作。

通常用事后估计方法选取步长 ,例如,记

为步长等于

的差商计算公式,给定误差

界 ,当

时,

就是合适的步长。

例 7.2 对函数

,取不同的步长计算

,观察误差变化规律,从而确定最佳步长。

解: 误差 0.10 0.09 0.08 0.07 0.06 3.1630 3.1622 3.1613 3.1607 3.1600 -0.0048 -0.0040 -0.0031 -0.0025 -0.0018 0.05 0.04 0.03 0.02 0.01 3.1590 3.1588 3.1583 3.1575 3.1550 误差 -0.0008 -0.0006 -0.0001 -0.0007 -0.0032

表中数据显示,当步长 从 0.10 减少到 0.03 时,数值微分误差的绝对值从 0.0048 减少到 0.0001,而 随着 的进一步减少,误差的绝对值又有所反弹,表明当步长 小于 0.03 时,舍入误差起了主要作用。

在实际计算中是无法得到误差的准确数值的,这时以 取 = 0.04。

最小为标准确定步长,本例中

7.1.2 插值型数值微分

对于给定的

的函数表,建立插值函数

,用插值函数

的导数近似函数

的导数。

设 多项式 ,以



上的节点,给定 的相应阶的导数,即

,以

为插值点构造插值

的各阶导数近似



时,

(7.5)

误差项为:

例 7.3 给定

,并有

,计算



解:作过

的插值多项式:



代入

得三点端点公式和三点中点公式:

利用泰勒(Taylor)展开进行比较和分析,可得三点公式的截断误差是



类似地,可得到五点中点公式和五点端点公式:

7.1.3 样条插值数值微分

把离散点按大小排列成 的样条函数 :

,用

关系式构造插值点







时,可用

计算导数。

7.2

数值积分

在微积分中用牛顿-莱布尼兹(Newton-Leibniz)公式计算连续函数

的定积分:

但是, 当被积函数是以点列

的形式给出时, 当被积函数

的原函数

难以得到时,例如 ,则无法用牛顿-莱布尼兹积分公式计算。有时当被积函数的原函数过于复 杂时,也不宜套用积分公式计算积分,而应采用数值积分公式计算定积分。 在微积分中,定积分是黎曼(Rimann)和的极限,它是分割小区间长度趋于零时的极限,即

在数值积分公式中,只能用有限项的和近似上面的极限,通常由函数在离散点函数值的线性组合形式 给出。

记 分值, 程。 称为积分节点,

,在本章中,用 称为积分系数。确定

表示精确积分值,用

表示近似积

中积分系数

的过程就是构造数值积分公式的过

怎样判断数值积分公式的效果?代数精度是衡量数值积分公式优劣的重要标准之一。

) 定义 7.1 (代数精度)记

上以

为积分节点的数值积分公式为

若 精度。

满足



, 则称

具有

阶代数

由此可知当

具有

阶代数精度时,对任意的

阶多项式都有



7.2.1 插值型数值积分

对给定的被积函数在 近似计算

上的点列 ,即

作拉格朗日插值多项式

,以



,则有

数值积分误差,也就是对插值误差的积分值



对一般的函数

,但若

是一个不高于 次的多项式,由于

,而有

。因此, 阶插值多项式型式的数值积分公式至少有 阶代数精度。

例 7.4 建立

上节点为

的数值积分公式。

解:由



得以数值积分公式:

7.2.2 牛顿-柯特斯(Newton-Cote’s)积分

把积分区间

分成 等分,记步长为 ,取

,取等分点

作为数值积分

节点,构造拉格朗日插值多项式

由此得到的数值积分称为牛顿—柯特斯积分。下面可以看到,牛顿-柯特斯积分系数和积分节点以及 积分区间无直接关系,系数固定而易于计算。 梯形积分





为插值节点构造线性函数

,有

那么,

提取公因子 没有任何关系了。

后,得到牛顿-柯特斯积分的组合系数:



,它们已与积分区间



(7.6)



为梯形积分公式。它的几何意义是用梯形面积近似代替积分值(图 7.2)。

图 7.2 梯形积分 怎样确定梯形积分公式的代数精度?

我们可以取

验证



时,有

即:



时,有

即:



时,有

得梯形求积公式具有一阶代数精度。 梯形积分公式的误差:





因为



上不变号,由积分中值定理得到梯形求积公式的截断误差:

(7.7)

Simpson 辛普森(Simpson Simpson)积分

对区间

作二等分,记 ,那么,有

。以



为插值节点构造二次插值函数

计算得到积分组合系数:







(7.8)

S (f)称为辛普森或抛物线积分公式。它的几何意义是用过三点的抛物线面积近似代替积分的曲边面积 (图 7.3)。

图 7.3 抛物线积分面积

分别将

代入到



中,可以得到

表明辛普森公式对于次数不超过三次的多项式准确成立, 一个三次多项式满足条件:

具有三阶代数精度。因此可设

, 计算得到误差为:

于是有



故辛普森求积公式的截断误差:

(7.9)

牛顿-柯特斯积分系数

等分区间

,取等分点为积分节点, 为插值节点构造插值函数 。

,其中

。以

其中



,代入上式得

(7.10)

这里称

为牛顿-柯特斯系数

可见在取等距节点时,积分系数

与积分节点和积分区间无直接关系,只与插值的节点总数有关,

而在例 7.3 中的积分系数是待定系数,这就简化了数值积分公式,而不必对每一组插值节点 xi 都要计算一 组相应的积分系数 ai。在公式(7.10)中取 =1,可算出梯形积分系数;取 =2,可算出辛普森积分系数。 在表 7.1 中列出 从 1 到 6 的牛顿-柯特斯系数。从表中可以看出牛顿-柯特斯系数具有对称性。 表 7.1

1 2 3 4

5 6

7.2.3 求积公式的收敛性与稳定性

定义 7.2 若

,则称求积公式

是收敛的。

定义中的

包含了

,通常都要求计算积分的求积公式是收敛的。

稳定性是研究计算公式 ,误差记为

当 。

有误差

时,

的误差是否增长,现设

定义 7.3 对任给

,只要

,就有

则称求积公式

是稳定的。

定理 7.1 若求积公式

的系数

,则求积公式是稳定的。

证明 由于

,

,故有

于是对

,只要

,就有

故求积公式是稳定的。

7.3

复化数值积分

由插值的龙格现象可知,高阶牛顿-柯特斯积分不能保证等距数值积分系列的收敛性,同时可证(略) 高阶牛顿-柯特斯积分的计算是不稳定的。因此,实际计算中常用低阶复化梯形等积分公式。

7.3.1 复化梯形积分

把积分区间分割成若干小区间,在每个小区间

上用梯形积分公式,再将这些小区间上的数值

积分累加起来,称为复化梯形公式。复化梯形公式用若干个小梯形面积逼近积分 比用一个大梯 形公式效果显然更好,如图 7.4 所示。这种作法使我们想起定积分定义,即它为被积函数无限分割的代数 和。这也正是计算定积分最朴素的算法。

图 7.4 复化梯形公式积分视图

复化梯形积分计算公式



作等距分割,有



,于是



上,

,则有

记 等分的复化梯形公式为

,有

(7.11)

复化梯形公式截断误差



,根据均值定理,当

时,存在

,有

,于是

(7.12)

由此看到复化梯形公式的截断误差按照

或者

的速度下降, 事实上, 可以证明, 只要



上有界并黎曼可积,当分点无限增多时,复化梯形公式收敛到积分





,则有

对于任给的误差控制小量

,有



就有

,式中

表示取其最大整数。

7.3.2 复化辛普森积分
把积分区间分成偶数等分 ,记 ,其中 是节点总数, 是积分子区间的总数。



, 。

,在每个区间

上用辛普森数值积分公式计算,则得

到复化辛普森公式,记为 复化辛普森积分计算公式



,称

(7.13)

为复化辛普森积分公式,它是



上采用辛普森积分公式叠加而得。下面用图 7.5 显示

复化辛普森积分计算公式中节点与系数的关系,取

,在每个积分区间上提出因子

后,三个节点的

系数分别是 1,4,1;将 4 个积分区间的系数按节点的位置累加,可以清楚地看到,首尾节点的系数是 1, 奇数点的系数是 4,偶数点的系数是 2。

1

4

1

1

4

1 1 4 1 1 4 4 1 1

1

4

2

4

2

4

2

图 7.5 复化辛普森积分系数

复化辛普森公式的截断误差



,在

上的误差为

因此,



(7.14)

与复化梯形公式类似,误差的截断误差按照

或者

的速度下降。可以证明,只要





有界并黎曼可积,当分点无限增多时,复化辛普森公式收敛到积分





,则有

对任给的误差控制小量

,只要



就有



例 7.5 求 多少?

, 计算中要求有 5 位有效数字。 用复化梯形和复化辛普森求积公式的分点应取

解:

由复化梯形误差公式得到:

计算出

,复化梯形公式至少要在

.00 等分 n = 68。

由复化辛普森误差公式,有

在复化辛普森公式中取





7.3.3 复化积分的自动控制误差算法
复化积分的误差公式表明,截断误差随分点 的增大而减小,对于给定的误差量 ,用估计函数导数 的界的方法可计算出 。用误差公式计算满足精度的分点数,像是在做一道计算导数 上界的微积

分习题 (如例 7.5 所示) 但是在实际运算中, 。 一般难以估计出函数的各阶导数界, 也就无法确定分点数 。 在计算中常用误差的事后估计方法,即用 估计误差 。

T2n (f )的计算公式

对定积分

,取分点

,计算得

取分点

,计算得

这里,

。可以看到,

的值是

与新增分点

的组合。

取分点

,计算得

这里,



同理, 计算

时只要在

的基础上计算新增分点



的值再做组合, 如图 7.6 所示。

图 7.6



一般地,每次总对前一次的小区间分半,分点加密一倍,并可充分利用老分点上的函数值,每次只需 计算新增分点的和。



上 等分,

,则有



上的中点为



(7.15)

其中





其中



类似地,可得积分节点为 ,

的辛普森求积公式的关系式:

(7.16)

其中:

由误差公式:

由于 ,于是

,分别为 及

个点上的均值,可视

上式表明

的误差大约是

误差的 4 倍。



(7.17)

由此得到启发, 对任给的误差控制量 而用

, 要

, 只需 在计算机上也不难实现。

即可,

作为控制手段简单直接,序列

复化积分的算法描述
从数值积分的误差公式可以看到,截断误差随分点 的增长而减少,控制计算的精度也就是确定分点 数 。在计算中不用数值积分的误差公式确定分点数 的理论模式,而用 通过增加分点自动满足精度的方法称为数值积分公式的自动积分法。即在计算中构造序列 作为控制, ,

直到



时停止计算,由分点数自动控制积分值的误差,并取



下面描述复化数值积分公式的自动控制误差算法,详细程序和算例请看本章 7.6 节。 1.输入:误差控制精度 e = eps;初始分点值 。

2.计算 分点的复化梯形积分

T1=T2+100

//迭代计算中 T1 和 T2 分别表示



3.while | T1-T2|>e T1 = T2

H=Hn

//计算新增节点的值

T2= (T1+H )/2 H = h/2,n =2n end while 4.输出积分值 T2。 在自动控制误差算法中初始分点值不宜过小,以防假收敛。 //将区间一分为二

7.3.4 龙贝格(Romberg)积分

由前面得到的关系式(7.17),可以将 即

作为

的修正值补充到



(7.18)

其结果是将梯形求积公式组合成辛普森求积公式,截断误差由

提高到

。这种手段称为外

推算法。外推算法在不增加计算量的前题下提高了误差的精度,是计算方法中一种常用手法。

不妨对

再做一次组合。由

得到

(7.19)

复化辛普森公式组成复化柯特斯公式,其截断误差是

。同理对柯特斯公式进行组合:

得到具有 7 次代数精度和截断误差是

的龙贝格公式:

还可以继续对

做上去。

为了便于在计算机上实现龙贝格算法,将 表示梯形、辛普森、柯特斯积分,行标 表示分点数 龙贝格计算公式:

统一用 或步长

表示,列标
j。

分别

对每一个

从 2 做到 ,一直做到

小于给定控制精度停止计算。

龙贝格算法

龙贝格算法按表 7.2 元素的行序进行运算,

在计算中每个元素只用到上一行和本行 ; 在

的元素。对上面的算法进一步优化,对每 k 行可将计算定义在两行元素之间,令 每计算一行元素后,要将 。 表 7.2 龙贝格算法计算元素顺序表

1.输入区间端点

,精度控制值 ,循环次数

,定义函数

,取



2.



3.for

=2 to

{

}

4.输出



7.4

重积分计算

在微积分中计算二重积分是用化为累次积分的方法进行的。计算二重数值积分也是计算累次数值积分 的过程。为了简化问题,我们仅讨论矩形域上的二重积分。有很多非矩形域上的二重积分可作变换将其转 换到矩形域上。

(7.20)

其中:

是常数,



上连续。像在微积分中一样,将二重积分化为累次积分:

(7.21)



二重积分的复化梯形公式

对区间



分别选取正整数

,在 轴和

轴上分别有步

用复经梯形公式计算

,计算中将 当作常数,有

(7.22)

再将

当作常数,在 方向上计算式(7.23)中每一项的积分,有

=



+

+

+

= 积分区域的 4 个角点的系数是 1/4,4 个边界的系数是 1/2,内部节点的系数是 1。 误差:



在积分区间内。

例子 7.6 用复化梯形公式计算二重积分

,取



解:

如表 7.3 所示:

表 7.3

数值表

的数值如表 7.4 所示:

表 7.4

数值表



=0.873601

的准确值是 0.886176。

二重复化辛普森求积公式

对区间



分别

等分和 等分,在 轴和

轴上分别有步长

均为偶数 类似于二重复化梯形公式推导,得



误差:



在积分区间内。

按例 7.6 的化分区间,

的值如表 7.5 所示:

表 7.5

7.5
对插值型积分公式

高斯(Gauss)型积分公式介绍

在牛顿-柯特斯积分公式中要求节点是等距的, 其优点是计算积分系数的公式规则相同, 缺点是制约了 求积公式的代数精度。可以证明:当节点个数 为偶数时,求积公式具有 为奇数时,求积公式具有 阶的代数精度。 阶的代数精度;当节点个数

如果我们不预先指定求积节点

的位置,

和权系数

都作为待定的常数,能否适当地确定它们, 个方程来确定,取一个函数组: 次的多项式,都可以用 次

以提高积分公式的代数精度?回答是肯定的。 ,这一组函数构成了

个待定常参数,需要

次多项式的基,任一小于等于

这组函数的线性组合来表示。如果某一积分公式,对这组函数都能精确积分,则此积分公式就有 代数精度。

例 7.7 计算求积系数 度。

和求积节点

, 使得

至少具有 3 阶代数精

解:按照求积公式的代数精度定义,分别令

,得方程组:

解方程组得:



求积公式:

按例 7.7 的方式,构造更高阶的代数精度的求积公式,生成求积系数和求积节点的方程组并无困难, 而求解该方程组则无一定的章法可循。一般地,通过计算正交多项式的零点作为求积节点。

当取积分节点为正交多项式的零点时,则节点个数是 的求积公式具有 分节点为正交多项式的零点的数值积分公式为高斯型积分公式。 为了一般性,考虑积分

阶的代数精度。并称积

其中

称为权函数。当

时,即是普通的积分。

对于不同的权函数

选定的节点也不相同。

如何构造高斯型积分公式呢?

对给定的

及权函数

,由施密特(Schmidt)正交化过程作出正交多项式 的 个零点 ,这个 个零点就是积分节点;以

;解出正交多项式 这些节点构造插值多项式,计算积分系数

其中

是拉格朗日插值基函数。

高斯型求积公式为

高斯型积分公式的优点是它的代数精度高,特别是对无穷区间或瑕积分更有效,但计算正交多项式的 零点即积分节点有一定的工作量,好在数值学家们已算出一些特定的函数的积分节点和积分系数,在计算 中我们可以查表直接得到这些数值。 本章并不构造各种高斯型积分公式,有关的详细内容请参考有关的教材。

下面给出

上,取权函数

的高斯型积分。

取[-1,1]上权函数

的正交多项式为勒让德(Legendre)多项式:

高斯-勒让德

相应的积分节点和积分系数表如下:

2

=-0.5773503, =-0.8611363,

=0.5773503 1.0000000,1.0000000 =-0.3399810 0.3478548,0.6521452 0.6521452,0.3478548

4 =0.3399810, -x1= x5=0.9061798 -x2= x4=0.538469 =0.8611363 0.2369269,0.2369269 0.4786287,0.4786287 0.5688889

5

=0.0

要计算一般区间

上的积分

,只需作变量代换

则有

,其中 求积,即

,这样

仍可用高斯积分

例 7.8 应用两点高斯-勒让德积分公式计算 解:



I = (-0.5773503)2cos (-0.5773503) + (0.5773503)2 cos (0.5773503)
=0.558608

例 7.9 应用两点高斯-勒让德积分公式计算 解:





,得到积分

例 7.10 证明不存在

使求积公式

的代数精度超过

次。

证明:只要能找到一个 节点 及

次多项式,使求积公式两边不相等即可。用反证法,假定存在求积系数和 使求积公式对任何 次多项式 精确成立。现取

代入求积公式左端得

而公式右端

,故右端与左端不相等,与假设矛盾,说明不存在

次代

数精度的求积公式,故高斯型求积公式是具有最高代数精度的求积公式。

7.6

程序示例

程序 7.1 复化梯形公式

的自动控制误差算法。

算法描述

输入

的值,定义





; T1=T2+100 while |T1-T2|>ε T1 = T2;

H=h

!或用复化梯形公式计算 T2

T2 = (T1 + H) /2;h= h/2;n = 2n; end while

输出 T2。 程序源码 /////////////////////////////////////////////// // Purpose:梯形公式的自动控制误差算法 ////////////////////////////////////////////// # include <stdio.h> # include <math.h> # define f (x) (sin (x)) # define m 100 # define a 1.0 # define b 2.0 # define epsilon 0.00001 void main( ) { int n = m; int i; double T n,H_n,T1,T2; double h = (b-a) / n; T_n = (f (a) + f(b)) /2; // f(x)由 define 定义 //

for (i =1;i<n;i+

T_n + = f(a+h

i);

T_n = h; T2 = T_n; T1 = T2+100;

do { T1 = T2;

for (i=0,H_n=0;i =1;i<n;i+

H_n+=f (a+h i+h/2);

H_n = h; T2 = (T1 + H_n) /2

h = h /2;n = n 2; } while (fabs (T1-T2)>epsilon) printf (“T=% 1f \ n”,T2); } // --------End of File---------

计算实例

对于

,在区间

上验证梯形公式的自动控制误差公式。

程序输入输出

对于不同

,修改程序# define f(x)项,本例

,初始

,区间

。)

T = 0.956447 程序 7.2 龙贝格积分算法 计算公式和算法描述第 7.3.4 节所述。

程序源码
////////////////////////////////////// // Purpose:龙贝格算法 //

////////////////////////////////////// # include<stdio.h> # include<math.h> # define f (x) # define N_H # define MAXREPT # define a # define b # define epsilon 1.0 2.0 0.00001 (sin (x)) 20 10

double compute T (double aa,double bb,long int n) // 复化梯形公式 { int i;double sum,h = (bb-aa)/n;sum = 0;

for (i=1;i<n;i+

sum + = f (aa + i h); sum + = (f (aa) + f (bb)) / 2;

return ( h sum); } void main ( ) { int i; long int n = N_H,m = 0; double T[2] [MAXREPT +1]; T[1][0] = compute T (a,b,n);

n = 2; for (m = 1;m<MAXREPT;m++) } for ( i = 0;i<m;i++ ) { T[0][i]=T[1][i];} T[1][0]=computeT(a,b,n);

n =2; for ( i =1;i<m + + =

T[1][i]=T[1][i-1]+(T[1][i-1]-T[0][i-1]) / (pow(2,2 m)-1); if ((T[0][m-1]-T[1][m])<epslon) { printf (“T=% lf \ n”,T[1][m]); return; } printf (“Return no solved… \ n”); } // ----------End of File--------- } // 输出

计算实例

利用龙贝格积分法计算



程序输入输出

(对于不同

,修改程序# define 。)

项,本例

,区间

,初始

对于本程序的给定,输出结果。 0.956447

第8章

常微分方程数值解

在描述系统的动态演变时,例如,物种的增长和蜕变,物体的运动,电路的振动瞬变,化学反应过程 等,都能将其表示为以时间 为变量的常微分方程或方程组。 物体冷却过程的数学模型可用下式表示:

它含有自变量 、未知函数 以及它的一阶导数

,是一个常微分方程。在微分方程中我们称只有一

个自变量函数的微分方程为常微分方程,自变量函数个数为两个或两个以上的微分方程为偏微分方程。给 定微分方程及其初始条件,称为初值问题;给定微分方程及其边界条件,称为边值问题。 本章主要讨论常微分方程的初值问题:

(8.1) 或记为

只有一些特殊形式的

,才能找到它的解析解;对于大多数常微分方程的初值问题,主要用数 在求解区间 表示 上剖分点列 的近似值。

值方法在计算机上求它的数值解。常微分方程初值问题的数值解是求 的数值解 。 在计算中约定

表示常微分方程准确解的值,

本章介绍求解微分方程的差分数值方法。解常微分方程初值问题的主要手段是差分方法,这是一种通用性 强适用面广的简单方法。

通常我们假定(8.1)中 有

满足李普希兹(Lipschitz)条件,即存在常数

,使对



则初值问题(8.1)的解存在且惟一。

8.1

欧拉(Euler)公式

8.1.1 基于差商的欧拉公式

对于常微分方程初值问题[式(8.1)],在求解区间

上作等距分割的剖分,步长

,记

。用数值微商的方法,即用差商近似微商数值求解常微分方程。

用向前差商近似

做出

的在

处的一阶向前差商式:



,于是得到



的近似值

可按

或 求得。类似地,由

以及

得到计算

近似值

的向前欧拉公式:

(8.2)

由差商(差分)得到的方程(8.2)称为差分方程。



直接算出

值的计算格式称为显式格式,向前欧拉公式是显式格式。

欧拉方法的几何意义

以 以

为斜率,通过点 为斜率过点

做一条直线,它与直线

的交点就是

。依此类推,



的直线与直线

的交点。欧拉法也称为欧拉折线法,如图 8.1 所示。

图 8.1 欧拉折线法 例 8.1 假定某公司的净资产因资产本身产生了利息而以 4%的年利率增长,同时,该公司以每年 100 万的数额支付职工工资。净资产的微分方程:

= 0.04

-100

(t 以年为单位)

分别以初始值 产趋势。

1500 万,

2500 万,

3500 万,用欧拉公式预测公司 24 年后的净资

解:



分别以

1500,

2500,

3500 代入,计算结果见表 8.1。 表 8.1 计算结果

1 1460

2500. 3540

5 1283.35 2500. 3716.65

2 1418.4 2500. 3581.6 6 1234.68 2500. 3765.32 3 1375.14 2500. 3624.86 7 1184.07 2500. 3815.93 4 1330.14 2500. 3669.86 8 1131.43 2500. 3868.57 9 1076.69 2500. 3923.31 17 552.1 2500. 4449.9

10 1019.76 2500. 3980.24 18 474.183 2500. 4525.82 11 960.546 2500. 4039.45 19 393.151 2500. 4606.85 12 898.968 2500. 4101.03 20 308.877 2500. 4691.12 13 834.926 2500. 4165.07 21 221.232 2500. 4778.77 14 768.324 2500. 4231.68 22 130.081 2500. 4869.92 15 699.056 2500. 4300.94 23 35.2845 2500. 4964.72 16 627.019 2500. 4372.98 24 63.3042 2500. 5063.3 从表 8.1 可以看到当利息赢利低于工资的支出,公司的净资产逐年减少,以至净资产为负值;当利息 赢利与工资的支出平衡时,公司的净资产每年保持不变;当利息赢利超过工资的支出,公司的净资产稳步 增长。在图 8.2 中 分别表示初始值 3500,2500 和 1500 的三条净资产趋势曲线。

图 8.2 三种初始值的净资产趋势

用向后差商近似

做出



处的一阶向后差商式:



,得到

的近似值

的计算公式:

类似地,可得到计算

近似值

的计算公式:

(8.3)

公式(8.3)称为向后欧拉公式。通常 程,称为隐式欧拉公式,需要通过迭代法求得



的非线性函数,因此式(8.3)是关于 。其中,初始值

的非线性方

可由向前欧拉公式提供。

从算法结构上看,显式公式比隐式公式简单;从方法的稳定性和精度上看,多数情况下隐式公式优于 显式公式。最简单的迭代公式(皮卡(Picard)迭代)为:

直到

给定精度。

可以证明,h 充分小时,以上迭代收敛。事实上,记

,则

充分小时,可以保证

,其中

为李普希兹(Lipschitz)常数。

用中心差商近似

做出

的在

处的中心差商式:



,可得到

的近似值

的计算公式:

类似地,可得到计算

近似值

的计算公式:

(8.4)

公式(8.4)称为中心公式。按公式(8.4),需要知道 用其它公式算出 采用。 ,再用中心格式算出

的值才能求得

的值。因此,要先

。需要指出的是,中心格式不是稳定格式,因此不能

8.1.2 欧拉公式的收敛性

定义 8.1 8.1(局部截断误差)在假设 称为局部截断误差。

,即

步是精确的前提下,考虑的截断误差





做泰勒展开:

(8.5)

欧拉公式(8.2)是由以上展开式中截断

而得,故欧拉公式的局部截断误差为:

(8.6)



为线性函数,则

,这时局部截断误差为 0,由欧拉公式得到的解为精确解,故欧拉

公式是一阶方法。

图 8.3 Euler 公式的误差 定义 8.2 如果给定方法的局部截断误差是

则称该方法是

阶的,或称具有

阶精度。

: 以下考查欧拉向后公式的局部截断误差:

故此方法的局部截断误差主项为

,也是一阶方法。

整体截断误差和收敛性

在计算 从计算

的局部截断误差时,是假定在



值是准确的前提下,即 的误差会扩散到

。实际上, 中,将这些

开始,每个

都会产生截断误差,而且在

前列点的误差累计到计算 这一误差。

的误差中,称为整体截断误差。它将影响到方法的收敛性,我们将估计

由微分方程理论,为保证微分方程解的存在惟一性及稳定性,通常 件,即对任意, ,存在 满足



应满足李普希兹条

于是由式(8.5)与(8.2)两式相减得到

曾记

,故有



,则有

因此有



,由公式

,最终可得:

其中

由原始误差引起,当初始值为精确值时,这一项的值是 0;

仍是由截断引

起,由于

,当

时,

,即欧拉公式是收敛的。

8.1.3 基于数值积分的差分方法

对常微分方程(8.1)即

两边在区间

上积分得:



用矩形近似积分公式计算



,有

则得向前欧拉公式:

(8.7)



,有

则得向后欧拉公式:

(8.8)

用梯形近似积分公式计算

则得到梯形公式:

(8.9)

梯形公式也是隐式格式,可用皮卡迭代或牛顿迭代计算(8.9)的

。计算中为了保证一定的精度,

又避免用迭代过程不菲的计算量,一般先用显式公式算出初始值,再用隐式公式进行一次修正。称为预估 -校正过程。例如,下面是用显式的欧拉公式和隐式的梯形公式给出的预估-校正公式:

(8.10) 式(8.10)也称为改进的欧拉公式,它可合并写成

例 8.2 分别用梯形公式和预估-校正公式解初值问题:

解:



用下面的迭代公式,对每个迭代 4 次,



该方程的精确是

,计算结果如表 8.2 所示。 表 8.2 计算结果

1 2 3 4

0.1 0.2 0.3 0.4

1.1118 1.2520 1.4331 1.6763

1.1111 1.2500 1.4236 1.6667

0.0007 0.0020 0.0095 0.0004

8.2

龙格-库塔方法

8.2.1 二阶龙格-库塔方法
常微分方程初值问题:



在 点的泰勒展开:

这里

。取

,就有

(8.11)

截断

可得到

近似值

的计算公式,即欧拉公式:

若取

,式(8.11)可写成:



(8.12)

截断

可得到

近似值

的计算公式:



上式为二阶方法,一般优于一阶的欧拉公式(8.2),但是在计算 点的值,因此,此法不可取。

时,需要计算



龙格-库塔设想用 (8.12)的主体,即用

在点



值的线性组合逼近式

(8.13) 逼近

得到数值公式:

(8.14) 或更一般地写成

对式(8.13)在

点泰勒展开得到:

将上式与式(8.12)比较,知当

满足

时有最好的逼近效果,此时式(8.13)-式(8.14) 有无数组解。

。这是 4 个未知数的 3 个方程,显然方程组

若取

,则有二阶龙格-库塔公式,也称为改进欧拉公式:

(8.15)

若取

,则得另一种形式的二阶龙格-库塔公式,也称中点公式:

(8.16)

从公式建立过程中可看到,二阶龙格-库塔公式的局部截断误差仍为

,是二阶精度的计算公式。 ,是四阶精

类似地,可建立高阶的龙格-库塔公式,同时可知四阶龙格-库塔公式的局部截断误差为 度的计算公式。

欧拉法是低精度的方法,适合于方程的解或其导数有间断的情况以及精度要求不高的情况,当解需要 高精度时,必须用高阶的龙格-库塔等方法。 四阶龙格-库塔方法应用面较广,具有自动起步和便于改变步长的优点,但计算量比一般方法略大。 为 了保证方法的收敛性,有时需要步长取得较小,因此,不适于解病态方程。

8.2.2 四阶龙格-库塔公式
下面列出常用的三阶、四阶龙格-库塔计算公式。

三阶龙格-库塔公式

(1)

(8.17)

(2)

(8.18)

(3)

(8.19)

四阶龙格-库塔公式

(1)

(8.20)

(2) 例 8.3 用四阶龙格-库塔公式(8.20)解初值问题:

(8.21)

解:取步长

,计算公式为:

计算结果列表 8.3 中。 表 8.3 计算结果

1

0.2

1.24789

1.24792

0.00003

2 3 4

0.4 0.6 0.8

1.63762 2.29618 3.53389

1.63778 2.29696 3.53802

0.00016 0.00078 0.00413

8.2.3 步长的自适应

欧拉方法和龙格-库塔方法在计算

时仅用到前一步

的值,我们称这样的方法为单步法。在单步 变

法计算中根据需要可以取等步长或变步长。 对于 化剧烈的区段,则步长可以取小一些。

变化平缓的区段, 步长可以取大一些; 而对于

以计算为

为例,已知

的数值,取

为初始步长, 是给定的控制精度,用欧拉或龙格-库塔

方法的计算公式。



为由

和取步长 算出的

,记

为由

和取步长

算出的 ,直到 ,直到

,即



均是

的近似值,当| 后取 为步长的值。而当

,逐步放大步长, ,逐步缩小步长,

时为止,最 时为止,

取 为步长的值。 下面给出自适应算法的描述:

给定控制精度 和初始步长



对于

IF

THEN

WHILE

ENDWHILE

ELSE

WHILE

ENDWHILE

ENDIF

8.3

线性多步法

常微分方程初值问题[式(8.1)]与积分 等价,在 8.1.3 节中用一些简单的数 值积分公式建立了欧拉和梯形等数值方法,下面给出一些更一般的用数值积分工具求解常微分方程初值问 题的方法。数值积分公式当积分节点中包含 点 ,有 时称为隐式公式,否则称为显式公式。分别取 为剖分

若用积分节点

构造插值多项式近似

,在区间

上计算数值积分

,则称为隐式公式。

与通常数值积分不同的是,在线性多步法公式中,有两个控制量 节点。

控制积分区间, 控制插值

特别取

,有

若积分

用关于积分节点

的数值积分近似, 就可得到显式的阿达姆斯公式。

例 8.4 构造

的显式公式:

这里



截断

可得到

近似值

的计算公式:

(8.22) 上式称为二阶显式阿达姆斯公式。类似可得 三阶显式阿达姆斯公式:

(8.23)

四阶显式阿达姆斯公式:

(8.24)

若积分 用关于积分节点, 式。隐式格式需要通过迭代求解。它们分别为 二阶隐式阿达姆斯公式(也称为梯形公式):

的数值积分近似, 就可得到隐式的阿达姆斯公

(8.25)

三阶隐式阿达姆斯公式:

(8.26)

四阶隐式阿达姆斯公式:

(8.27)

例 8.5 建立 解:

显式公式。

计算格式:

其中: 截断误差:



这是一个三步三阶的显式格式,可用龙格-库塔公式计算



例 8.6 构造

显式公式。

解: 取

为积分区间, 以

构造拉格朗日插值多项式。

[



计算格式: 截断误差:

为了避免迭代,可用预估-校正公式替换上面的计算公式。

8.4

常微分方程组的数值解法

8.4.1 一阶常微分方程组的数值解法
将由 个一阶方程组成的常微分方程初值问题:

(8.28) 写成向量形式:

(8.29) 其中:







前面对常微分方程所用的欧拉方法、 龙格-库塔方法等各种方法, 都可以平行地应用到常微分方程组的 数值解中。为了叙述方便,下面以两个方程组为例,给出相应的计算公式。 常微分方程组:

欧拉公式:

(8.30)

预估-校正公式:

(8.31) 四阶龙格-库塔公式:

(8.32)

例 8.7 两种果树寄生虫,其数量分别是

,其中一种寄生虫以吃另一种寄生虫为生,

两种寄生虫的增长函数如下列常微分方程所示,预测 3 年后这一对寄生虫的数量。

解:记

在本题中

。用欧拉预估-校正公式(8.31):



,计算结果见表 8.4。 表 8.4 计算结果 (年) 1 2 3 4 1.6 0.28766 0.291308 0.00021727 1.2 1.33103 1.47776 1.63947

8.4.2 高阶常微分方程数值方法
下面以三阶常微分方程为例说明高阶常微分方程的数值计算步骤。

(8.33)

将三阶方程化为一阶方程。令

(8.34) 得到一阶方程组:

已将高阶方程化为一阶方程组了

8.5

常微分方程的稳定性

用欧拉公式

计算

时,假设计算中的某一步有误差,而以后的计算是精确的,

那么这一步的误差对以后的计算有何影响? 若某算法在计算过程中任何一步产生的误差在以后的计算中都逐步衰减,称该方法是绝对稳定的;如 果这一步的误差在以后的计算中恶性放大,则称方法是不稳定的。 在选用方法时,方法的绝对稳定性是最为重要的指标。不绝对稳定的方法是不能采用的。讨论绝对稳 定性是把方法用到最为典型的稳定微分方程(8.35)上进行的。

(8.35) 例 8.8 对方程(8.35)讨论欧拉方法的稳定性。

解:用欧拉方法计算方程(8.35)的公式:



有误差

,记为

,那么

就有误差

,记为



满足 将(8.36)与欧拉公式相减得到:

(8.36)



(8.37)



时,即

落在图 8.4 中的单位圆上时有

即这时误差逐次衰减,方法是绝对稳定的。 图 8.4 上的单位圆,称为绝对稳定区域。

图 8.4 欧拉方法的绝对稳定区域 例 8.9 对方程(8.35)讨论向后欧拉方法的稳定性。 解:用向后欧拉方法计算方程(8.35)的公式:

误差方程:

计算相邻两步误差的比值:



, 恒有

, 当绝对稳定区域是左半平面时, 则称在这个区域数值方法是

稳定的,或称无条件绝对稳定。 例 8.10 对方程(8.35)讨论中心差分方法的稳定性。 解: 用中心差分方法计算方程(8.35)的公式:

误差方程:

(8.38)

不失一般性,我们考虑由



对以后

的影响,差分方程(8.36)的特征方程为:

它的两个根为:



由差分方程理论(略),

的一般解可由下式表达:

(8.39)

其中





决定。

由于

,故

因 增大而恶性发展,所以方法为不稳定的。 可以证明,龙格-库塔方法以及显式或隐式的阿达姆斯方法均是绝对稳定的。

8.6

程序示例

程序 8.1 用改进的欧拉公式来求解常微分方程初值问题:

算法描述

对给定的

,用改进的欧拉公式

求解常微分方程初值问题的解。

程序源码
////////////////////////////////////

// purpose:改进的区拉公式 /////////////////////////////////// # include <stdio.h> # include <math.h>

//

# it main ( ) { int m; int i; double a,b,y0; double xn,yn,xnl,ynl,ynlb; double h; printf (“nInput the begin and end of x:”); scanf (“%lf%lf ”,&a,&b); printf (“Input the y value at %f:”,a); scaf (“%lf ”,&y0 ); printf (“nInput m value [divide (%f,%f):”,a,b]; scanf(”%lf ”,&m); if ( m<= 0) { printf (“Please input a number larger then 1. \ n”); return 1; } h = (b-a) / m; xn = a;y = y0; for (i =1;i<= m;i++)

{ xnl = xn + h;

ynlb = yn + h f (xn,yn);

ynl = yn + h /2 (f (xn,yn) + f (xnl,ynlb)); printf (“x %d = %f,y %d = %f \ n”,i,xnl,i,ynl); xn = xnl;yn = ynl; } return 0;}

计算实例
用改进的欧拉公式,求解常微分方程初值问题的解:

程序输入输出
Input the begin and end of x: 0 0.4 Input the value at 0.000000:1 Input m value [divide (0.000000,0.400000)]:4 x1=0.100000,y1=1.110500 x2=0.200000,y2=1.248276 x3=0.300000,y3=1.424760 x4=0.400000,y4=1.658736 程序 8.2 用四阶龙格-库塔方法求解常微分方程初值问题:

算法描述

对给定的

,用四阶龙格-库塔方法求解常微分方程初值问题。

程序源码
//////////////////////////////////////////////////////// // purpose:四阶龙格-库塔法求解常微分方程初值问题 // //////////////////////////////////////////////////////// # include <stdio.h> # include <math.h>

int main ( ) { int m; int i; double a,b,y0; double xn,yn,ynl; double k1,k2,k3,k4; double h; printf (“\ nInput the begin and end of x:”); scanf (“%lf %lf ”,&a,&b); printf (“Input the y value at %f:”,a);

scaf (“%lf ”,&y0 ); printf (“nInput m value [divide (%f,%f):”,a,b]; scanf(”%d”,&m); if ( m<= 0) { printf (“Please input a number larger then 1. \ n”); return 1; } h = (b-a) / m; xn = a;yn = y0; for (i =1;i<= m;i++) { kl=f (xn,yn);

k2=f ( (xy + h /2),(yn + h kl /2));

k3=f ( (xn + h / 2),(yn + h k2 /2);

k4=f ( (xn + h),(yn + h k3);

yn1=yn + h/6 (k1+2 xn + = h;

k2+2 k3+k4);

printf (“x % d =% f,y %d = % f \ n”,i,xn,i,ynl); yn = ynl; } return 0; }

计算实例
用四阶龙格-库塔公式解初值问题:

程序输入输出
Input the begin and end of x: 2.0 2.6 Input the value at 0.000000:1 Input m value [divide (2.000000,2.600000)]:3 x1=2.200000,y1=1.356505 x2=2.400000,y2=1.661361 x3=2.600000,y3=1.939104



推荐相关:

数值计算方法1

数值计算方法1 - 天津理工大学计算机科学与技术专业数值计算方法实验一报告。仅供参考,请勿抄袭~


数值计算方法复习提纲

数值计算方法复习提纲_理学_高等教育_教育专区。西南交通大学数值计算方法复习提纲 数值计算方法复习提纲 第一章 数值计算中的误差分析 1.了解误差及其主要来源,误差...


数值计算方法教案

数值计算方法教案 - 《计算方法》教案 课程名称: 适用专业: 适用年级: 任课教师: 编写时间: 计算方法 医学信息技术 二年级 张利萍 2011 年 8 月 新疆医科大学...


数值计算方法期末复习答案终结版

数值计算方法期末复习答案终结版_理学_高等教育_教育专区 暂无评价|0人阅读|0次下载|举报文档数值计算方法期末复习答案终结版_理学_高等教育_教育专区。一、 名词...


《数值计算方法》试题集及答案

数值计算方法》试题集及答案_财会/金融考试_资格考试/认证_教育专区。《数值计算方法》复习试题 一、填空题: ? ? ? 4 ?1 0 ? ? A ? ? A?? ? 1 4...


数值计算方法

数值计算方法_工学_高等教育_教育专区。1.已知 ln(2.0)=0.6931;ln(2.2)=0.7885,ln(2.3)=0.8329, 试用线性插值和抛物插值计算.ln2.1 的值并估计误差...


数值计算方法第一章答案

数值计算方法第一章答案 - 数值计算方法第一章答案 2.(1)3580 绝对误差限: ;相对误差限: 经过四舍五入得到的近似值 2580,其各位都是有效数字,故有四 位...


数值计算方法-简单迭代

数值计算方法-简单迭代 - 《数值计算方法》实验 2 报告 班级: 姓名: 学号: 成绩: 1. 实验名称 实验 2 非线性方程的迭代解法(之简单迭代法) 2. 实验题目 ...


数值计算方法

数值计算方法_数学_自然科学_专业资料。本科实验报告 课程名称: 计算机数值方法 实验项目: (名称) 实验地点: 行勉楼 B208 专业班级: 软件 1401 学号: 学生姓名...


数值计算方法

数值计算方法_数学_自然科学_专业资料。\要求: 1. 独立完成,作答时要写明题型、题号; 2. 作答方式:手写作答或电脑录入,使用 A4 格式白纸; 3. 提交方式:...

网站首页 | 网站地图
All rights reserved Powered by 学霸学习网 www.tceic.com
copyright ©right 2010-2021。
文档资料库内容来自网络,如有侵犯请联系客服。zhit325@126.com