当前位置:首页 >> 高中教育 >> 第二章 一元非线性方程的解法

第二章 一元非线性方程的解法


计算方法 吴筑筑编

一元非线性方程的解法
孙剑 计算机学院信息管理系
1

本章主要内容:
?2.1

初始近似根的确定

?2.2
?2.3 ?2.4 ? ?

二分法
迭代法的一般知识 牛顿迭代法(切线法) 弦截法(割线法)

埃特金(Aitken)迭代法
2

2.5 2.6

引言:
? 本章讨论:求一元非线性方程f(x)=0的根的近似 值(其中f(x)为非线性函数) 。
? 例如:x8-x3+5x-3=0 ——f(x)为多项式,则

f(x)= 0为代数方程。
? 例如:x-sin

x=0, e-x-x=0 ——f(x)含有三角函数、 指数函数或其他超越函数,则f(x)= 0为超越方 程。

3

引言:
? 如果

? 本章讨论:求一元非线性方程f(x)=0的根的近似值

f(x*)=0 ,则x*是方程f(x)=0的根,也称它 是函数f(x)的零点。

? x*为方程f(x)=0的m重根(函数f(x)的m重零点)

当且仅当f(x*)=f ′(x*)=f″(x*)=…=f(m-1)(x*)=0, 但 f(m)(x*)≠0成立,m是正整数。
? 方程的

1 重根称为单根,这时 f(x*)=0 而

f′(x*)≠0。
4

引言:
? 求解代数方程: x8-x3+5x-3=0,

超越方程:x-sin x=0, e-x-x=0的精确解存在一定 的困难。 ? 求根的大致步骤:
? (1)判定根的存在性。 ? (2)确定根的初始近似值(初始近似根)。 ? (3)根的精确化。

5

2.1 初始近似根的确定

? 定理 如果函数f(x)在区间[a,b]上连续,严格单调, 且f(a)f(b)<0,则在(a,b)内方程有且仅有一个实根。
y f(x) b x 0 a x* ? 假设f(x)在某区间(a,b)内有且仅有一个实根x*,若 b – a较小,则可在(a,b)上任取一点x0作为初始近 似根。一般情形,可用逐步扫描法等。
6

2.1 初始近似根的确定
? 例2-1:f(x)=x3-x-1=0 由于 f(0)= -1<0, f(∞)>0,
? 故方程至少有一正实根。设从x=0 出发,取

h=0.5为步长向右扫描,得可见在(1,1.5)内 有根。 x 0 0.5 1.0 1.5 f(x) ― ― ― +
2 ? f ( x) ? 3x ?1 ? 0 (1 ? x ? 1.5)

? 所以f(x)在区间[1,1.5]上单调连续,因而在

(1,1.5)内有且仅有一个实根,故可在[1 ,1.5]上 任取一点做初始近似根。 7

2.1 初始近似根的确定
? 求初始近似根的逐步扫描法:(设(a,b)内有根) (1)x0←a; (2)若f(x0)f(x0+h)≤0则结束;否则做下步。 (3)x0← x0+h,转(2) (其中h为预选的步长)
x0 x0x +h x +h x +h x +h 0 x0 0 x0 0 x0 0 x0+h a

8

2.2 二分法
? 2.1.1
? 2.1.2

二分法的基本思想和计算过程
近似根的误差估计

? 2.1.3
? 2.1.4

二分法的计算流程
二分法举例

9

2.1.1 二分法的基本思想和计算过程
? 条件:函数f(x)在[a,b]上连续,严格单调,且 f(a)f(b)<0,则方程f(x)=0在区间(a,b)内有且 仅有一个实根x*。 ? 二分法的基本思想
? 将含方程根的区间平分为两个小区间,然后判

断根在哪个小区间,舍去无根的区间,而把有 根的区间再一分为二,再判断根属于哪个更小 的区间,如此周而复始 ,直到求出满足精度 要求的近似根。
10

2.1.1 二分法的基本思想和计算过程

x*

11

2.1.1 二分法的基本思想和计算过程
? 具体计算过程:
1 第1次二分,取中点 x0 ? (a ? b) 2 若 f(a)f(x0 )<0,则 x*∈( a, x0 ), 否则x*∈(x0 ,b )

令a1=a , b1=x0或令a1=x0 , b1=b新的有根区间为
(a1 , b1 ) ,长度是原来的一半。
y a1 0 a

f ( x)

y
0 a

f ( x) a1 x0 b1 x* b x
12

x*

b1 x0

b x

2.1.1 二分法的基本思想和计算过程
1 ? 第2次二分,取中点 x1 ? ( a1 ? b1 ) 2 ? 若 f(a1 )f(x1 )<0,则 x*∈( a1 , x1 ),令a2=a1 , b2=x1;否则 令 a2=x1 , b2=b1新的有根区间为 (a2 , b2 ) 。
y
0 a

f ( x) a1 x0 a2
x1 b2

b1 b x

13

2.1.1 二分法的基本思想和计算过程
? 如此反复,有(a, b) ? (a0 , b0 ) ? (a1, b1) ? (a2 , b2 ) ? ? (ak , bk ) ?

1 xk ? (ak ? bk ) ? (ak , bk ),k=0,1,2 2
1 1 bk ? ak ? (bk ?1 ? ak ?1 ) ? k (b ? a ) 2 2

?0 (k ? ?)

ak ? xk ? bk

lim ak ? lim xk ? lim bk ? x
k ?? k ?? k ??

*

14

2.1.1 二分法的基本思想和计算过程
? 近似根xk的误差估计:

1 | x ? xk |? (bk ? ak ) 2
?

? bk ?1 ? ak ?1
1 ? k ?1 (b ? a ) 2

x* l m |—————|—————| ak xk bk

15

2.1.1 二分法的基本思想和计算过程
? 由此得二分过程的结束原则:
? 先给定精度要求ε(绝对误差限)
? (1)当|bk+1 – ak+1|<



ε时结束二分计算,取 k+1 ,取

16

x*≈xk ;
? (2)事先由ε估计出二分的最小次数

x*≈xk
*

b?a 由 x ? xk ? k ?1 ? ? 2 ln(b ? a) ? ln ?
k ?1 ? ln 2

得2

k ?1

?

b?a

?

2.1.1 二分法

? 例2-2:求方程f(x)=x3-x-1=0在区间(1, 1.5)内的根。 要求用二分法,精确到10 -2。 ? 解: a=1, b=1.5, f(1)= -1<0 , f(1.5)>0.
*

0.5 b?a ?2 由 x ? xk ? k ?1 ? ? , 有 k ?1 ? 0.5 ? 10 , 2 2

2

k ?1

? 10

2

? 得k+1 =7 满足要求, 二分7次计算到 x6 即可。 1 (1 ? 1.5) f(1.25)<0, f(1)f(1.25)>0, 2 得有根区间[ 1.25,1.5];

得有根区间[ 1.25,1.375];

1 (1.25 ? 1.5) ? 1.375, 2

f(1.375)>0, f(1.25)f(1.375)<0,
17

2.1.1 二分法的基本思想和计算过程
如此继续,计算结果见列表: x6=1.3243, 根 x*≈1.32
k+1 ak bk xk f(xk)的符号

1
2 3

1
1.25 1.25

1.5
1.5 1.375

1.25
1.375 1.312 5


+ -

4
5 6 7

1.312 5 1.37 5
1.312 5 1.343 8 1.312 5 1.328 2 1.320 4 1.328 2

1.343 8
1.328 2 1.320 4 1.324 3


+ - -
18

注意:
?―精确到10p‖ 绝对误差限ε=0.5×10 p
绝对误差限 ε=10 p |bk+1 – ak+1|< ε

?―精度要求为ε=10p‖

19

2.1.1 二分法的基本思想和计算过程
? 二分法 的计算 流程:
开始

输入a bε f(a)
(a+b)/2

? 实验: P13例2-1

y
x

x

b



yf(x)<0 b-a<ε 是 输出x 结束



x

a


20

练习:
?P25 习题1

21

2.3 迭代法的一般知识
?2.3.1 迭代法的基本思想及几何意义
?2.3.2 迭代法的收敛条件及误差估计式 ?2.2.3 局部收敛性与迭代法收敛的阶

22

2.3.1 迭代法的基本思想及几何意义

?1.

迭代法的基本思想

?2.

迭代法的几何意义

23

2.3.1 迭代法的基本思想:
? 方程f(x)=0化为等价形式的方程 x=g(x)

? 构造迭代公式 xk+1= g(xk ) , k=0,1,2,……
? 取初始近似根x0 ,迭代计算

x1=g(x0),x =g(x ),…….. 得到迭代序列 {x }, x 2 1 k 例 f(x)=xe -1=0,可化为等价形式
–xlim ? 当g(x)(称为迭代函数)连续, x=e –x g(x)= 若

e

lim xk ?1 ? lim g ( xk ) k ??
k ??

迭代公式 xk+1=e –xk
*

k ??

xk ? x* ,则由

x ? g(x )
*
24

等价地有f(x*)=0, 故x* 即为方程的根。

2.3.1 迭代法的基本思想:
? 实际计算到 |xk – xk-1 |<ε(ε是预定的精度),取 x*≈xk 。 ? 迭代法也称逐次逼近法。 ? 迭代序列 {xk}= x1 ,x2 ,x3 ,……xk…… ? 迭代公式收敛(发散)指迭代序列 {xk} 收敛(发 散)。 问题:迭代公式是否一定收敛?

25

2.3.1 迭代法的基本思想:
? 例2-3:求方程 f(x)=x – 10x+2=0的一个根,取4位 有效数字计算。 ? 解:因 f(0)= 1 > 0, f(1)= -7< 0, 所以方程在(0,1)中有根。方程改写为两种等 价形式:

(1) x ? g1 ( x ) ? lg(x ? 2)
对应的迭代公式分别为:

(2)

x ? g2 ( x) ? 10x ? 2
xk

(1) xk ?1 ? lg( xk ? 2); (2) xk ?1 ? 10 ? 2.

26

2.3.1 迭代法的基本思想:
? 用迭代公式(1) 取 x0=1, 算得
xk ?1 ? lg( xk ? 2)

x1=lg3=0.4771, x2=lg(x1+2)=0.3939,

…., x6=0.3758, x7=lg(x6+2)=0.3758,… x6、x7重合 所以迭代公式(1)是收敛的,x*≈0.3758。

27

2.3.1 迭代法的基本思想:
? 用迭代公式(2) xk ?1 ? 10xk ? 2 x1=10-2=8, x2=108-2≈108,
8 10 x3=10 -2≈ 8 10 10 ,……

迭代公式(2)发散。

28

2.3.1 迭代法的几何意义:
迭 代 法 的 几 何 意 义

迭代函数g(x)满足什么条件时,迭代序列才收敛? 29

2.3.2 迭代法的收敛条件 及误差估计式

30

2.3.2 迭代法的收敛条件及误差估计式
? 定理2-1 设迭代函数g(x) 满足: (1)任意x∈[a, b]时 a≤ g(x) ≤b (2) g(x)可微( g(x) 在[a, b]上有连续的一阶导 数),存在正数 q < 1,使对任意 x ∈[a, b] 都有 |g′(x)| ≤ q < 1 则方程x=g(x)在 [a, b]上有惟一的根x* ;

且对于[a,b]上任意初始近似根x0 ,迭代公式 xk+1=g(xk)均收敛于方程的根x*;还有误差估计式
k q q * x ? xk ? xk ? xk ?1 ? x1 ? x0 1? q 1? q
31

2.3.2 迭代法的收敛条件及误差估计式
? 证明: 先证根的存在性。 记 f(x)=x-g(x), 由条件(1)有 a≤g(a)、g(b)≤b , 于是 f(a)=a-g(a)≤0, f(b)=b-g(b)≥0,

由于f(x)是连续函数,故必存在 x* ∈[a, b] 使
f(x*)=0 于是x*为方程 x=g(x)的根。

32

2.3.2 迭代法的收敛条件及误差估计式
其次,证根的惟一性 : 设y*也是方程的根,则x*=g(x*), y*=g(y*), |x*- y*|=|g(x*) – g(y*)| = |g′(ξ)(x*- y*)| 微分中值定理, ξ在x*,y*之间 ≤q |x*- y*|

由已知条件 |g ′(x)| ≤q <1

故必有x*=y*,所以方程在[a,b]的根是惟一的。
33

2.3.2 迭代法的收敛条件及误差估计式
再证迭代公式收敛 由x 0∈[a, b] 知 xk∈[a, b],k =0,1,2,…… |x*-xk | =|g(x*)-g(xk-1) |= |g′(ξ)(x*-xk-1)|
k→0 0<q<1, q ≤q | x*-xk-1 | =q|g(x*)-g(xk-2) |= q|g′(η)(x*-xk-2)|

≤q2|x*-xk-2| ≤q3| x*-xk-3|≤……

lim xk ? x
k ??

*

≤qk |x*-x0| →0 (k→∞)

所以,对[a,b]上任取的x0 ,迭代公式xk+1= g(xk )都收 敛于x*。 34

2.3.2 迭代法的收敛条件及误差估计式
? 定理2-1 设迭代函数g(x) 满足: (1)任意x∈[a, b]时 a≤ g(x) ≤b (2) g(x)可微( g(x) 在[a, b]上有连续的一阶导 数),存在正数 q < 1,使对任意 x ∈[a, b] 都有 |g′(x)| ≤ q < 1 则方程x=g(x)在 [a, b]上有惟一的根x* ;

且对于[a,b]上任意初始近似根x0 ,迭代公式 xk+1=g(xk)均收敛于方程的根x*;还有误差估计式
k q q * x ? xk ? xk ? xk ?1 ? x1 ? x0 1? q 1? q
35

2.3.2 迭代法的收敛条件及误差估计式

36

2.3.2 迭代法的收敛条件及误差估计式

37

2.3.2 迭代法的收敛条件及误差估计式
? 定理2-2’ 设在区间[a, b]上方程 x=g(x) 有根x*,且 对一切x∈[a, b] 都有|g′(x)| ≥ 1,则对于该区间 上任意x0(≠x*), 迭代公式xk+1= g(xk )一定发散。 证明

xk ?[a, b] 且x0 ? x 时,
*
* ? xk ? x ?| g (? ) | xk ?1 ? x *

? xk ?1 ? x* ? ...... ? x0 ? x* ? 0

xk ? x

*

不可能收敛于0。
38

2.3.2 迭代法的收敛条件及误差估计式

39

2.3.2 迭代法的收敛条件及误差估计式

40

2.3.2 迭代法的收敛条件及误差估计式
? 判断迭代公式收敛,需要迭代函数g(x) 满足: (1)任意x∈[a, b]时 a≤ g(x) ≤b (2) g(x)可微( g(x) 在[a, b]上有连续的一阶导 数),使对任意 x ∈[a, b] 都有 |g′(x)| < 1

? 判断迭代公式发散,需要判断迭代函数g(x) :
|g′(x)| ≥1
41

2.3.2 迭代法的收敛条件及误差估计式
? 定义2-1 若存在x*的某个邻域D:|x-x*|<δ,使得 则称xk+1=g(xk)在根x*附近具有局部收敛性。

迭代公式xk+1=g(xk)对于任意的初值x0∈D均收敛,

42

2.3.2 迭代法的收敛条件及误差估计式
? 定理2-2 设在区间 [ a , b ]上方程 x = g(x)有根 x*,且g(x)有连续的一阶导数。如果有正数 q<1 , 使得对任意 x∈[a, b] 都有|g′(x)| ≤ q < 1,则存在 x*的某个邻域,只要x0 属于此邻域,迭代公式 xk+1= g(xk )必收敛于x*。(也称迭代公式有局部 收敛性) x*-δ x*+δ * x a b
43

2.3.2 迭代法的收敛条件及误差估计式
? 证明:

存在x*的邻域(x*-δ, x*+δ)?[a, b],
当x∈ (x*-δ, x*+δ)时,

|g′(x)| ≤ q < 1,

| x-x* | <δ,

由中值定理得|g(x)-x*|=|g(x)-g(x*)|

=|g′(ξ)(x-x*)|≤q| x-x* |<δ,
于是g(x) ∈ (x*-δ, x*+δ),这说明在该邻域内定理 的条件(1)也成立。所以结论成立。
44

2.3.2 迭代法的计算步骤
迭代法的计算步骤: (1)确定方程 f(x)=0 的一个等价形式 x=g(x), 要

求在某个含根区间[a,b]内满足① a≤ g(x) ≤ b

②|g′(x)| ≤ q < 1 。
(2)构造迭代公式 xk+1=g(xk),选取初始近似根x0 ,
进行迭代计算。

(3) 当|xk – xk-1 |<ε(ε是预定的精度)时停止计算,

取x*≈xk 。
45

2.3.2 迭代法的收敛条件及误差估计式
开始 迭代 法的 计算 框图 输入x0ε g(x0) x1

|x1-x0 | <ε 是 输出x1 结束



x1

x0

46

2.3.2 迭代法的收敛条件及误差估计式
–x在x=0.5附近的一个根,按5位 求方程 x=e 例2-4: 小数计算,结果的精度要求为ε=10 –3。
–x=0. 方程等价于 f(x)=x – e 解: 由于f(0.5)<0, f(0.7)>0, 故x*∈(0.5, 0.7), 令g(x)=e –x,g'(x)= -e –x, 0.5≤x ≤0.7时, 0.5≤ g(x) ≤0.7

在(0.5, 0.7)内, g(x) 的一阶导数连续,且有

| g?( x) |?| (e )? |? e ? e
?x ?x

?0.5

? 0.61 ? 1

所以用迭代公式 xk+1=e –xk 进行计算是收敛的。 47

2.3.2 迭代法的收敛条件及误差估计式
x0=0.5, x1=e –x0=0.60653, x2=e –x1=0.54524,……. 迭代结果:
k
0 1 2 3 4 5

xk+1=e –xk

xk
0. 5 0. 606 53 0. 545 24 0. 579 70 0. 560 07 0. 571 17

xk – xk-1
0. 106 53 -0. 061 29 0. 034 46 -0. 019 63 0. 011 10

k
6 7 8 9 10

xk
0. 564 86 0. 568 44 0. 566 41 0. 567 56 0. 566 91

xk – xk-1
-0. 006 31 0. 003 58 -0. 002 03 0. 001 15 -0. 000 65

|x10 - x9 |=0.00065<ε, 故 x*≈x10≈ 0.567

48

2.3.2 迭代法的收敛条件及误差估计式
3-2x-5=0在(1.5, 2.5)内有一实根, 方程 f(x)=x 例2-5 分别取

作为迭代函数,试判断对应的迭代公式是否 收敛,并用一个收敛的迭代公式求方程的根, 精度要求为ε=10 -4。
49

1 3 (2) x ? g 2 ( x) ? ( x ? 5) 2

(1) x ? g1 ( x) ? 3 2x ? 5

2.3.2 迭代法的收敛条件及误差估计式
? 解: (1) 当1.5 ? x ? 2.5时, g1 ( x ) 单调,且

g1 (1.5) ? 2, g1 (2.5) ? 3 10 ? 2.15443.....

故 1.5 ? g1 ( x) ? 2.5,
2 ? ( x) ? ? g1 3 1 (2 x ? 5)
2 3

2 ? ? 3

1 (2 ?1.5 ? 5)
2 3

2 1 1 ? ? ? 3 4 6

所以迭代公式 收敛。

xk ?1 ? g1 ( xk ) ? 3 2 xk ? 5
50

2.3.2 迭代法的收敛条件及误差估计式
取x0=2 计算,
结果列表:
k 0 1 2 3 4 5

x1 ? 3 2 x0 ? 5 ? 3 9 ? 2.08008

x2 ? 3 2 x1 ? 5 ? 2.09235
xk 2 2.08008 2.09235 2.0942 2.09449 2.09454 xk-xk-1 0.08008 0.01227 0.00187 0.00027 0.00005

因为 |x5-x4|=0.00005<10-4 所以x*≈x5=2.0945

51

2.3.2 迭代法的收敛条件及误差估计式
1 3 (2)迭代函数 g 2 ( x ) ? ( x ? 5) 2
3 2 2 ? | g 2 ( x) |? x ? 1.5 ?1.5 ? 1 x ? 2.5时, 2

当1.5 ?

根据定理2-2, 迭代函数对应的迭代公式

1 3 xk ?1 ? ( xk ? 5) 2
是发散的。
52

2.2.3 局部收敛性与迭代法收敛的阶
p阶收敛:若xk→x*并且存在p≥1,c>0,使

lim

xk ?1 ? x * xk ? x *
p

k ??

?c

线性收敛p=1(0<c<1)

超线性收敛2>p>1
平方收敛p=2

53

2.2.3 局部收敛性与迭代法收敛的阶
定理2-4:若g ( x)在x ? g ( x)的根x *附近有连续的p ( p ? 1)阶 导数,且g '(x* )=g"(x* )= ? g(p-1)(x* )=0,g(p)(x* ) ? 0。 则迭代公式xk+1 ? g ( xk )在x *的邻域是p阶收敛的,

且 lim

xn ?1 ? x * xn ? x *
p

n ??

?

g

( p)

( x*) ?0 p!

54

2.4 牛顿迭代法(切线法)
?2.4.1
?2.4.2 ?2.4.3

牛顿迭代法的基本思想
牛顿迭代法的收敛性 牛顿迭代法的计算步骤

55

2.4 牛顿迭代法的基本思想
? 牛顿迭代法的基本思想:
? 把非线性方程线性化,用线性方程的解逐步逼

近非线性方程的解。

56

2.4 牛顿迭代法的基本思想
①过曲线上的点pk(xk , f(xk))作切线,取切线与x轴

的交点为 xk+1。

②切线的方程 y=f(xk)+f ?(xk)(x – xk) ,点(xk+1 , 0)满 足该方程,即 0= f(xk)+f ?(xk)(xk+1 – xk)
f ?(xk)(xk+1 – xk) = - f(xk)

若 f ?(xk )≠0,则得Newton迭代公式 f ( xk ) xk ?1 ? xk ? ,k ? 0,1,2,? f ?( xk )

57

2.4 牛顿迭代法的基本思想

x3 x2

6

58

2.4 牛顿迭代法的收敛性
? 牛顿法迭代函数:

f ( x) g ( x) ? x ? f ?( x)

f '( x) f '( x) ? f ( x) f ??( x) | f ??( x) | | g ?( x) |?|1 ? |? | f ( x) | 2 2 [ f ?( x)] [ f ?( x)]

?若x*是单根(即f(x*)=0, f ?(x*)≠0), f ″(x)连续从而有界,则当x足够靠近 x*时, 就能满足 |g′(x)|≤q <1, 保证迭代公式局部收敛。
59

2.4 牛顿迭代法的基本思想
? 定理2.3 如果在有根区间 [a,b]上 f′(x)≠0,f″(x)

f ( x0 ) f ??( x) ? 0 则牛顿迭代法收敛。

连续且不变号,在 [a,b]上取初始近似根 x0 ,使得

60

2.4 牛顿迭代法的计算步骤
(1)给出x0 , ε;

(2)计算

x1 ? x0 ?

f ( x0 ) f ?( x0 )

(3)若 x1 ? x0 ? ? 则转(4);否则 x0 ? x1 ,转(2); (4)输出x1 ,结束。

61

2.4 牛顿迭代法的计算步骤
例:用牛顿迭代法求方程 xex-1=0 在x=0.5附近的 根(取5位小数计算),精度要求为ε=10 –3。 解:

f ( x) ? xe ?1
x

f ?( x) ? e ? xe
x
xk

x

xk e ? 1 相应的牛顿迭代公式为 xk ?1 ? xk ? xk xk e ? xk e ? xk xk ? e ? xk ? 取x0=0.5,经计算可得: 1 ? xk

x ? x3 ? 0.567

?

62

2.4 牛顿迭代法的基本思想
? 例:用牛顿迭代法计算 3 解:令 x ? 3 ,则x2-3=0,求 3 等价于求方程 f(x)=x2-3=0的正实根。因为 f′(x)=2x ,由牛顿迭 2 代公式得 xk ?3 1 3
xk ?1 ? xk ? 2 xk ? ( xk ? ) 2 xk k ? 0,1, 2,

取初值 x0=1.5,迭代5次可得 3 ≈1.732050808 问题:如何用牛顿法计算任意正数的算术平方根? 是否还能用牛顿法计算一个正数的立方根?
63

2.4 Newton法变形牛顿
? 简化Newton法.
? 为减少计算导数的花费,可只求f ′(x0)以后所

有导数不另求。这相当于第一次作切线,以后 作其平行线。
f ( xk ) xk ?1 ? xk ? , k ? 0,1, 2, f ?( x0 )
?当然,这样收敛要慢些。还可以取折衷方案,

隔几步计算一下导数。
64

2.4 Newton法变形牛顿
? Newton-下山法.
? 每次迭代在改变量前加一因子以保证收敛:
?

f ( xk ) xk ?1 ? xk ? ?n , k ? 0,1, 2, f ?( xk )

?

这儿λn在0~1间,可用各种方法搜索,例如用分 半法取1,1/2,1/4,…试探,使 ∣f(xk+1)∣<∣f(xk)∣
65

2.4 牛顿迭代法的基本思想
牛顿迭代公式:
f ( xk ) xk ?1 ? xk ? ,k ? 0,1,2,? f ?( xk )
? 定理2.3 如果在有根区间 [a,b]上 f′(x)≠0,f″(x) 连续且不变号,在 [a,b]上取初始近似根 x0 ,使得
?( x) ? 0 f ( x0 ) f ? 则牛顿迭代法收敛。

66

2.5 弦截法(割线法)
? 研究目的:在牛顿法基础上,构造既有较高的收 敛速度,又不须导数的迭代公式。 f ( xk ) ? f ( xk ?1 ) ? 思想: 用差商 代替导数 f ?( xk ) xk ? xk ?1 弦截迭代公式
xk ?1 f ( xk ) ? xk ? ( xk ? xk ?1 ),k ? 1,2,? f ( x k ) ? f ( x k ?1 )
67

2.5 弦截法(割线法)

68

2.5 弦截法(割线法)
用弦截迭代法求上一节的方程 xex-1=0 在x=0.5附 近的根。 解方程化为 x-e –x=0, 令 f ( x) ? x ? e? x

弦截迭代公式为
*

x ? x3 ? 0. 567

xk ? e? xk xk ?1 ? xk ? ( xk ? xk ?1 ) ? xk ? xk ?1 ( xk ? xk ?1 ) ? (e ? e ) xk k
0 1 2 3 0.5 0.6 0.567 54 0.567 15
69

开始
输入x0 , x1 ,?

弦截法 的计算 框图:

f ( x1 ) x 2 ? x1 ? ( x1 ? x0 ) f ( x1 ) ? f ( x0 )
| x 2 ? x1 |? ? F x0 ? x1 x1 ? x 2

T
出 x 出 输 2

结束

70

2.6 埃特金(Aitken)迭代法
研究目的:有时迭代过程收敛缓慢,如何加速? 埃特金迭代法的构造条件:

已有迭代公式 xk+1= g(xk) 假定g′(x) 在求根区间内 改变不大 。
埃特金迭代公式

? 迭代 ? ? 改进 ? ?

~ ? g( x ),x ? g( x ~ ) x k ?1 k k ?1 k ?1 2 ~ ( x k ?1 ? xk ?1 ) xk ?1 ? xk ?1 ? ~ ?x xk ?1 ? 2 x k ?1 k

71

2.6 埃特金(Aitken)迭代法
例用埃特金迭代法求 x3-x-1=0 在 (1, 1.5)内的根。
解将方程改写成 x=x3- 1,建立迭代公式

xk ?1 ? x ?1, k ? 0,1, 2,
3 k

埃特金迭代公式为? x ? x 3 ? 1, x ? x 3 ? 1 k ?1 k k ?1 k ?1

? ( xk ?1 ? xk ?1 ) 2 ? ? xk ?1 ? xk ?1 ? x ? 2 x ? x k ?1 k ?1 k ?

x0 ? 1.5
72

2.6 埃特金(Aitken)迭代法
x0 ? 1.5
k
0 1 2 3 4 5

xk
2.375 00 1.840 92 1.491 40 1.347 10 1.325 18

xk
12.396 5 5.238 88 2.317 28 1.444 35 1.327 14

xk
1.5 1.416 29 1.355 65 1.328 95 1.324 80 1.324 72

73

74

计算方法 吴筑筑编

第三章 线性代数方程组的解法
孙剑
计算机学院信息管理系
75

本章主要内容:
? 3.1 ? 3.2 ? 3.3 ? 3.4 ? 3.5 ? 3.6 ? 3.7 ? 3.8 ? 3.9

简单迭代法的一般形式

雅可比迭代法和高斯—赛德尔迭代法
超松弛(SOR)迭代法 顺序高斯消去法 选主元高斯消去法 行列式和逆矩阵的计算 追赶法 三角分解法 线性方程组的最小二乘解
76

引言:
? 本章解决的问题:线性代数方程组的基本解法。 n个未知量n个方程的线性代数方程组

(3-1)
矩阵形式: AX=b, 其中:

当A非奇异时,det A≠0,方程组有唯一解。

77

引言:
两类解法: 迭代解法:构造适当的向量序列,用某种极限过 程去逐步逼近精确解。

例如:雅可比迭代法、高斯赛德尔迭代法等。
直接解法:经过有限步运算就能求得精确解的方 法,但实际计算中由于舍入误差的影响,这类方 法也只能求得近似解。 例如:高斯消去法、三角分解法等。
78

引言:
? 线性代数方程组的迭代解法 : 1、 简单迭代法的一般形式 2 、雅可比(Jacobi)迭代法

3 、高斯-赛德尔(Gauss-Seidel)迭代法 4 、超松弛(SOR)迭代法

79

3.1 简单迭代法的一般形式
? 线性代数方程组的迭代解法的基本思想:
? 构造适当的迭代公式,任选一个初始向量x(0)

进行迭代计算,使生成的向量序列。

x ,x ,

(0)

(1)

x

(k )

收敛于方程组的精确解x(*)。

80

3.1 简单迭代法的一般形式
方程组 Ax=b , A非奇异, 化为等价的方程组。

x ? Mx ? g
? m11 ?m M ? ? 21 ? ? ? mn1 m12 m22 mn 2 m1n ? m2 n ? ? ? ? mnn ?

(3-2)
? g1 ? ?g ? g ? ? 2? ? ? ? ? ? gn ? ? x1 ? ?x ? x ? ? 2? ? ? ? ? ? xn ?

,

,

81

3.1 简单迭代法的一般形式
方程组 Ax=b , A非奇异, 化为等价的方程组。

x ? Mx ? g
? x1 ? ? m11 ? x ? ?m ? 2 ? ? ? 21 ? ? ? ? ? ? ? xn ? ? mn1 m12 m22 mn 2 m1n ? ? m2 n ? ? ? ? mnn ?

(3-2)
? x1 ? ? g1 ? ?x ? ?g ? ? 2? ? ? 2? ? ? ? ? ? ? ? ? ? xn ? ? g n ?
82

3.1 简单迭代法的一般形式
? x1 ? m11 x1 ? m12 x2 ? ? m1n xn ? g1 ? ? x2 ? m21 x1 ? m22 x2 ? ? m2n xn ? g 2 ? ? ?x ? m x ? m x ? ? m x ? g n1 1 n2 2 nn n n ? n

xi =

?

n

mij x j + gi , i = 1, 2,L , n,
83

j= 1

3.1 简单迭代法的一般形式
构造简单迭代公式

x
x
(0)

( k ? 1)

? Mx

(k )

? g , k ? 0,1,2,?

(k ) (k ) (k ) (k ) T x = [ x , x , L , x (k = 0,1,L ), 其中 1 2 n ]

任取。M为迭代矩阵.
( k ?1 ) ? x1 ? ( k ?1) ? x2 ? ? ? x ( k ?1) ? n

? m11 x

(k ) 1

? m12 x

(k ) 2

? ? ? m1 n x

(k ) n

? g1

(k ) (k ) (k ) ? m 21 x1 ? m 22 x 2 ? ? ? m2n xn ? g2

????
(k ) (k ) (k ) ? m n 1 x1 ? mn2 x2 ? ? ? m nn x n ? gn

xi( k + 1) =

?

n

mij x(jk ) + gi , i = 1, 2,L , n,

(k = 0,1, 2, L )
84

j= 1

3.1 简单迭代法的一般形式
若存在极限
k

lim xi( k ) = xi*
(k )

,

i = 1, 2,L , n

* * * * T { x } 则称向量序列 收敛于向量 x = [ x1 , x2 ,L , xn ] (k ) * lim x = x 并记为

k

这时,称简单迭代法

x

( k + 1)

= Mx

(k )

+g

,

k = 0,1, 2,L

是收敛的,否则就是发散的。
85

3.1 简单迭代法的一般形式
例3-1 用迭代法求解线性方程组 ì 2 x + x = 3 ? 1 2 ? í ? ? ? - 2 x1 + 5 x2 = 3 解 构造方程组的等价方程组
ì x1 = - x1 - x2 + 3 ? ? í ? ? ? x2 = 2 x1 - 4 x2 + 3

据此建立迭代公式
(0) =0 取 x1(0) = x2 计算得 :

( k + 1) (k ) (k ) ì ? x = x x 1 1 2 + 3 ? í ( k + 1) (k ) (k ) ? x = 2 x 4 x 1 2 + 3 ? ? 2

(3) (4) (5) 祆 镲 x = 9 x = 15 x 1 1 1 = 33 镲 L 眄(3) (2) (4) (5) x2 = - 3, 镲 x2 = 9, x2 = - 15, x2 = 33, 镲 铑 迭代解离精确解 x1 = 1, x2 = 1 越来越远迭代不收敛 86

(1) 祆 镲 x 1 = 3 镲 眄(1) 镲 x2 = 3, 镲 铑

x1(2) = - 3

3.1 简单迭代法的一般形式
? 定义 设有矩阵序列{A(k)},其中A(k)=(aij(k))n,
k=1,2,.…, 若A(k)的每个元素序列{aij(k)}分别收敛 于A*的对应元素aij , 则称 {A(k)}收敛于矩阵A* ? 定理 矩阵序列{A(k)}收敛于矩阵A*的充要条件 是存在一种矩阵范数||· ||使 lim A( k ) ? A* ? 0
k ??

87

3.1 简单迭代法的一般形式
? 定理3-1 简单迭代公式 x(k + 1) = Mx( k ) + g , k = 0,1, 2,L

收敛的充要条件是迭代矩阵M的谱半径 r (M ) < 1
证:必要性 设迭代公式收敛,当k→∞时,

x

(k )

? x

*

则在迭代公式两端同时取极限得 x* = Mx* + g
k

lim x
(k )

( k + 1)

= lim Mx
k

(k )

+ g

记e

(k )

=x - x则 e

*

(k )

收敛于0 (零向量)
88

3.1 简单迭代法的一般形式
记 e( k ) = x( k ) - x* 则 e ( k ) 收敛于0 (零向量)
e(k ) = x(k ) - x* = Mx( k- 1) + g - (Mx* + g ) = M ( x( k- 1) - x* ) = Me( k- 1)

于是 由于 e

e
(0)

(k )

= Me

( k - 1)

= M e

2 ( k - 2)

=L = M e

Mk - 0
Mk ? 0

k (0)

可以是任意向量,故 e

(k )

收敛于0当且仅
0

k M 当 收敛于零矩阵,即当 k

矩阵序列:M1,M2,M3……Mk 收敛于零矩阵

89

3.1 简单迭代法的一般形式
于是 0 ? (r (M )) 所以必有
k
k . M r (M )

k

0

r (M ) < 1

充分性: 设 r (M ) < 1 , 则必存在正数ε, 使 M ? r (M ) + e < 1
0# M
k

k

M
k

k

? (r ( M )

e)

k

k

lim (r ( M ) + e ) = 0

k

所以 lim M

= 0即 k
*

lim M k = 0

故 e( k ) = x( k ) - x* 收敛于0
90

x ( k )收敛于 x

3.1 简单迭代法的一般形式
? 定理3-1 简单迭代公式 x(k + 1) = Mx( k ) + g , k = 0,1, 2,L 收敛的充要条件是迭代矩阵M的谱半径 r (M ) < 1

91

3.1 简单迭代法的一般形式
?定理3-2(迭代法收敛的充分条件) 若迭代矩阵M的一种范数 M < 1 ,则迭代公式
x( k + 1) = Mx( k ) + g
M 1- M
k

(k = 0,1,L )

收敛,且有误差估计式,且有误差估计式
x - x
* (k )

?

x( k )

x( k- 1)



x - x

*

(k )

M ? 1- M

x (1)

x (0)

92

3.1 简单迭代法的一般形式
收敛时令k→∞,有 等价地有Ax*=b . 控制迭代结束的实用标准:

x = Mx + g

*

*

max x
1? i? n

(k ) i

?x

( k ?1) i

??

其中?是指定的精度要求。 解向量

x* ? x

(k )
93

3.1 简单迭代法的一般形式
0.5 x2 ? 4 ? x1 ? ?0.2 x1 ? x2 ? 3 例3-1 ? 改写为 ? x ? ?0.2 x ? 3 ? 2 1 ? x1 ? 0.5 x2 ? 4
(0) ? 0 ,构造迭代公式 取 x1(0) ? x2

( k ?1) (k ) ? x ? 0.5 x ? 1 2 ?4 ? ( k ?1) (k ) x ? ? 0.2 x ?3 ? ? 2 1

k ? 0,1, 2,.....
, k ? 0,1, 2,

其矩阵形式为 x( k ?1) ? Mx( k ) ? g 其中

0.5? ? 0 ? 4? M?? , g?? ? , ? ? ?0.2 0 ? ? 3?

x

(k )

?x ? ?? ? ?x ?
(k ) 1 (k ) 2
94

3.1 简单迭代法的一般形式
( k ?1 ) (k ) ? x ? 0 . 5 x ?4 1 2 根据迭代公式 ? ? ( k ?1) (k ) ? ? ?0.2 x1 ?3 ? x2

计算得

(1) ? x ? 1 ?4 ? (1) ? ? x2 ? 3

(2) (1) ? x ? 0.5 x ? 1 2 ? 4 ? 0.5 ? 3 ? 4 ? 5.5 ? (2) (1) x ? ? 0.2 x ? 3 ? ?0.2 ? 4 ? 3 ? 2.2 ? ? 2 1

? ? x ? 0.5 x ? 4 ? 0.5 ? 2.2 ? 4 ? 5.1 , ? (3) (2) ? 3 ? ?0.2 ? 5.5 ? 3 ? 1.9 ? ? x2 ? ?0.2 x1
(3) 1 (2) 2

设精度要求为ε=0.005

,计算结果:
95

3.1 简单迭代法的一般形式
k 0
0

1
4

2

3

4
4.95

5

6

7

x

(k ) 1

5.5 5.1

4.99 5.005 5.001

x

(k ) 2

0

3

2.2 1.9

1.98

2.01 2.002 1.999
(7) 2

由于

x

(7) 1

?x

(6) 1

?? , x

?x
(7) 2

(6) 2

??

故求得方程组的解为

x1 ? x

(7) 1

? 5.001 , x2 ? x

? 1.999
96

?3.2 雅可比迭代法 和高斯—赛德尔迭代法

97

3.2.1 雅可比迭代法
?公式的构造:方程组Ax=b, 设A非奇异,其主对

角元aii≠0,(i=1,2,…,n),并且绝对值相对来说比较
大。 从第i个方程 ai1x1 ? ai 2 x2 ?
? aii xi ? ? ain xn ? bi

中分离出xi,得到等价的方程组

n 1 xi ? (bi ? ? aij x j ),i ? 1,2,? , n aii j ?1 j?i

雅可比迭代公式:
xi
( k ?1) j?i

n 1 (k ) ? (bi ? ? aij x j ), aii j ?1

( i ? 1,2,? , n)
98

3.2.1 雅可比迭代法
? x3 ? 9 ?10 x1 ? ??2 x1 ? 10 x2 ? x3 ? 7 精度要求为?=0.005。 ? ? x2 ? 5 x3 ? 4 ?

例3-2 求解方程组

解:等价的方程组为:

构造雅可比迭代公式

? x3 ) /10 ? x1 ? (9 ? ? x3 ) /10 ? x2 ? (7 ? 2 x1 ? x ? (4 ? x2 )/5 ? 3

(k ) ? x1( k ?1) ? (9 ? x3 ) /10 ? ( k ?1) (k ) (k ) x ? (7 ? 2 x ? x ? 2 1 3 ) /10 , ? ( k ?1) (k ) x ? (4 ? x )/5 2 ? 3

k ? 0,1, 2,
99

3.2.1 雅可比迭代法
迭代计算

x(0) ? 0 ? [0,0,0]T

( k ?1) (k ) ? x1 ? (9 ? x3 ) / 10 ? ( k ?1) (k ) (k ) x ? ( 7 ? 2 x ? x ? 2 1 3 ) / 10 ? ( k ?1) (k ) x ? ( 4 ? x 2 )/ 5 ? 3

x ? 0.9
(1) 1

x ? 0.7
(1) 2

x ? 0.8
(1) 3

(1) ? x1(2) ? (9 ? x3 ) /10 ? 0.98 ? (2) (1) (1) ? x2 ? (7 ? 2 x1 ? x3 ) /10 ? 0.96 ? (2) (1) x ? (4 ? x ) / 5 ? 0.94 3 2 ?

计算结果:
100

3.2.1 雅可比迭代法
k 0 1 2 3 0.99 4 5 …

x1(k) 0 0.9 0.98 0.994 0.9992 0.99980
x2(k) 0 0.7 0.96 0.998 0.99964




x3(k) 0 0.8 0.94 0.992
(5) (4) x2 ? x2 ? 0.00164 ? 0.005

0.998 0.99960



x1(5) ? x1(4) ? 0.0006 ? 0.005

x

(5) 3

?x

(4) 3

? 0.0016 ? 0.005

所以x=x(5)满足精度要求.
101

3.2.1 雅可比迭代法
问:上例中雅可比迭代公式的迭代矩阵是什么?
(k ) (k ) ? x1( k ?1) ? (9 ? x3 ) /10 ? 0.1x3 ? 0.9 ? ( k ?1) (k ) (k ) (k ) (k) x ? (7 ? 2 x ? x ) /10 ? 0.2 x ? 0.1 x ? 2 1 3 1 3 ? 0.7 ? ( k ?1) (k ) (k ) x ? (4 ? x ) / 5 ? 0.2 x ? 0.8 2 2 ? 3

答:由雅可比迭代公式

可以看出迭代 矩阵是:

? 0 ? B ? ?0 .2 ? ? 0

0 0 0 .2

0 .1? ? 0 .1? 0 ? ?

102

3.2.1 雅可比迭代法
xi

( k ?1)

雅可比迭代公式的分量形式:

1 (k ) ? (bi ? ? aij x j ), ( i ? 1,2,?, n) aii j ?1
j?i

n

xi

( k ?1)

a i i ?1 ( k ) a i i ? 1 ( k ) ai n ( k ) bi ai 1 ( k ) ai 2 ( k ) ? ? x1 ? x2 ? ....? x i ?1 ? xi ?1 ? ....? xn ? aii aii aii aii aii aii
i=1,2,…,n
103

3.2.1 雅可比迭代法
雅可比迭代的矩阵形式为

x
? ? 0 ? ? ? a21 ? a22 B?? ? ? ? ? ? a n1 ?? ? ann

( k ?1)

? Bx
a13 ? a11 a23 ? a22 ? an 3 ? ann

(k )

?f
a1n ? ? a11 ? ? a2 n ? ? a22 ? , ? ? ? ? ? ? 0 ? ? ? b1 ? ?a ? ? 11 ? ? b2 ? ? a22 ? f ?? ? ? ? ? ? ? ? ? bn ? ?a ? ? nn ?
104

a12 ? a11 0

? ? ?

an 2 ? ann

?

B为雅可比迭代矩阵。

练习:
?1 2 ?2? ? ? (1) 系数矩阵A ? ?1 1 1 ? , ? ?2 2 1 ? ?

? 0 ?2 2 ? ? ? 则B ? ? ?1 0 ?1? , ? ? ?2 ?2 0 ? ?

105

3.2.2 高斯-赛德尔迭代法
( k ?1) (k ) ? x1 ?思想:设想把 x1(k+1) 算出后立即代替 x1(k)用于后面 ? (9 ? x3 ) / 10

? ( k ?1) (k ) (k ) (k+1) x ? ( 7 ? 2 x ? x ? 2 x2 分量的计算,当 算出后立即代替 x2(k)用于后面 1 3 ) / 10 ? ( k ?1) (k ) x ? ( 4 ? x 3 ,期望这样会收敛得快些。 2 )/ 5 ?… 分量的计算,

? x3 ? 9 ?10 x1 例3-3 用高斯-赛德尔 ? ? ?2 x1 ? 10 x2 ? x3 ? 7 迭代法求解方程组 ? ? x2 ? 5 x3 ? 4 ? ? x3 ) /10 ? x1 ? (9 方程组化为等 ? ? x3 ) /10 ? x2 ? (7 ? 2 x1 价的方程组 ? x ? (4 ? x2 )/5 ? 3 106

3.2.2 高斯-赛德尔迭代法
k ?1) ( k ?(1) (( kk)) 高斯 -赛德尔迭代公式 雅可比迭代公式: ? x ? (9 ? x ) /10 /10 ? x1 1 ? (9 ? x33 ) ( k ?1) ( k ?1) (k ) ? ? ( k ? 1) ( k ) ( k ) ) /10 x ? (7 ? 2 x ? x ? 迭代计算: 2 1 3 ? (7 ? 2 x1 ? x3 ) /10 ? x? 2 ( k ?1) ( k ?1) (0) T x ? (4 ? x )/5 ? (k ? 1) (k 3 2) ? x ? 0 ? [0,0,0] x3 ? (4 ? x2 )/5 ? (1) x1 ? 0.9 (1) x2 ? (7+2×0.9+0)/10 = 0.88 (1) x3 ? (4+0.88)/5 = 0.976

x

(2) 1

(1) x2

? (9+0.976)/10 = 0.9976 ? (7+2×0.9976+0.976)/10 = 0.99712
? (4+0.99712)/5 = 0.99942
107

x

(1) 3

3.2.2 高斯-赛德尔迭代法
计算结果: k 0 1 0.9 2 0.9976 0.99712 3 0.99994 0.99993 0.99999 … … … … x1(k) 0

x2(k) 0 0.88

x3(k) 0 0.976 0.99942

x

(3) i

?x

(2) i

? 0.005 ,

i ? 1, 2,3

故x(3)即为满足精度的近似解,得 x1=0.999 94 , x2=0.999 93 , x3=0.999 99
108

3.2.2 高斯-赛德尔迭代法

高斯-赛德尔(Seidel)迭代公式 i ?1 n 1 ( k ?1) ( k ? 1) (k ) xi ? (bi ? ? aij x j ? ? aij x j ), (i ? 1,2,?, n) aii j ?1 j ? i ?1 aij bi , fi ? , i, j ? 1, 2, , n (i ? j ) 令 bij ? ? aii aii 迭代公式即
?x ? 0 ? b x ? ? b x ? f1 ? ( k ?1) ( k ?1) (k ) ? x2 ? b21 x1 ? 0 ? ? b2 n xn ? f 2 ? ? ? x ( k ?1) ? b x ( k ?1) ? b x ( k ?1) ? ? 0 ? f n1 1 n2 2 n ? n
( k ?1) 1 (k ) 12 2 (k ) 1n n
109

3.2.2 高斯-赛德尔迭代法
其矩阵形式是

x
bij ? ? aij aii

( k ?1)

? Lx
bi fi ? , aii

( k ?1)

? Ux

(k )

?f



i, j ? 1, 2,

, n (i ? j )

?0 ?b 21 ? L? ? ? ? ?bn1

0 bn 2

? ?0 b12 ? b1n ? ? ? ? 0 ? b 2n ? ? ,U ? ? ? ? ? ? ? ? ? ? ? ? 0? 0? ?
110

雅可比迭代矩阵B与L、U的关系为: B=L+U

3.2.1 雅可比迭代法
雅可比迭代公式:
xi
( k ?1) n 1 (k ) ? (bi ? ? aij x j ), aii j ?1 j?i

( i ? 1,2,? , n)

? ? 0 ? ? ? a21 ? a22 B?? ? ? ? ? ? a n1 ?? ? ann

a12 ? a11 0

a13 ? a11 a23 ? a22 ? an 3 ? ann

? ? ?

an 2 ? ann

?

a1n ? ? a11 ? ? a2 n ? ? a22 ? , ? ? ? ? ? ? 0 ? ?

? b1 ? ?a ? ? 11 ? ? b2 ? ? a22 ? f ?? ? ? ? ? ? ? ? ? bn ? ?a ? ? nn ?

111

3.2.2 高斯-赛德尔迭代法
(i ? 1,2,?, n)

高斯-赛德尔(Seidel)迭代公式 i ?1 n 1 ( k ?1) ( k ? 1) (k ) xi ? (bi ? ? aij x j ? ? aij x j ), aii j ?1 j ? i ?1 其矩阵形式是

x

( k ?1)

? Lx

( k ?1)

? Ux

(k )

?f

?0 ?b L ? ? 21 ? ? ? ?bn1

0 bn 2

? ?0 b12 ? b1n ? ? ? ? 0 ? b2 n ? ? ,U ? ? ? ? ? ? ? ? ? ? ? ? 0? 0? ?

112

3.2.2 高斯-赛德尔迭代法
改写迭代公式为: x(k+1)- L x(k+1)= U x(k)+f ,



(I– L) x(k+1)=U x(k)+ f , I:单位矩阵

x(k+1)=(I – L ) -1U x(k)+ (I – L ) –1f 可见,这也是简单迭代法,其迭代矩阵为

M = (I – L ) -1U
(高斯-赛德尔迭代矩阵)
113

3.2.2 高斯-赛德尔迭代法
前面的例题中,高斯-赛德尔迭代公式为:
?x ? (9 ? ( k ?1) ( k ?1) x ? (7 ? 2 x ? 2 1 ? ( k ?1) ( k ?1) x ? (4 ? x 2 ? 3
( k ?1) 1

? x ) /10 ? ( k ?1) (k ) ? 0.2 x ? x3 ) /10 1
(k ) 3

(k ) 0.1x3 ? 0.9

(k ) ? 0.1x3 ? 0.7

)/5 ?

( k ?1) 0.2 x2

? 0.8

高斯-赛德尔迭代矩阵为:
?0 ? ?1 0 0 ? ? 0 0 0? ? ?0 0 0.1? ? ?? ? ? ? ? ? ? M ? ? ?0 1 0? ? ?0.2 0 0? ? ?0 0 0.1? ? ? 0 ? ?0 0 1 ? ? 0 0.2 0? ? ?0 0 0 ? ? ? ? ?? ? ? ?0 ??
?1

0.1 ? ? 0 0.12 ? 0 0.024 ? ? 0
114

开 始

输 入 A, b ,y,? , N

0

,k ? 0

雅可 比迭 代法 的计 算框 图:
F

k ? k ?1 x ? y

k ? N0
F
yi

T

? (? ? a i j x
j ?1 j?i

n

j

? bi ) / a ii

i ? 1, 2 , ? n

1?i ? n

max yi ? xi ? ?
T

输 出 y1 , y 2 , ? , y n

输出失败信息

结 束

115

开 始

输入A, b, x , ? , N 0 , k ? 0

高斯赛得 尔迭 代法 的计 算框 图

k ? k ?1, d ? 0

T

k ? N0
F

对 i ? 1 , 2 , ? , n 计算 s ? x
i

x

i

?

(?

?

n

a

j ?1 j ? i

ij

x

j

? bi ) ? s }

a

ii

d ?

max{

d , x

i

d ??
T
输出失败信息

F

输出 x 1 , x 2 , ? , x n
结 束

116

3.2.2 高斯-赛德尔迭代法
GS迭代法和Jacobi迭代法的比较:

通常当两种方法都收敛时,GS迭代法往往收敛快 一些。但有时Jacobi迭代法收敛而GS迭代法发散; 有时又相反。

117

3.2.2 高斯-赛德尔迭代法
例如方程组
?0.2 x1 ? 0.5 x2 ? 2 ? ?0.4 x1 ? 0.5 x2 ? 3
( k ?1) (k ) ? x ? ? 2.5 x ? 1 2 ? 10 ? ( k ?1) (k ) x ? ? 0.8 x ?6 ? ? 2 1

雅可比迭代公式

发散,

但GS迭代公式收敛。 方程组
?4 x1 ? 5 x2 ? 30 ? ?2 x1 ? 5 x2 ? 20

的雅可比迭代公式
( k ?1) (k ) ? ? 1.25 x2 ? 7.5 ? x1 ? ? ( k ?1) (k ) x ? ? 0.4 x ?4 ? ? 2 1

收敛而GS迭代公式发散。

118

练习:
?1 2 ?2? ? ? (1) 系数矩阵A ? ?1 1 1 ? , ? ?2 2 1 ? ?

? 0 ?2 2 ? ?, 则B ? ? ? 1 0 ? 1 ? ? ? ? ?2 ?2 0 ? ?

所以

? 0 0 0? ?0 ?2 2 ? ? , U ? ?0 0 ?1? L?? ? 1 0 0 ? ? ? ? ? ? ? ?2 ?2 0 ? ? ?0 0 0 ? ?

? ?I ? B ? 1
2

2

?2 1 ? ?3

?
2

?

因而J迭代 法收敛。
119

? ( B) ? 0 ? 1

练习:
?而
? I ? B1 ? ? I ? ( I ? L) ?1U

? ? ( I ? L)?1 ( I ? L) ? ( I ? L) ?1U ? ( I ? L) ?1 ? ( I ? L) ? U

? 2 ?2 ? ? ( I ? L) ? U ? ? ? 1 ? ? (? ? 2) 2 2? 2? ? ?1 ? 0, ?2 ? ?3 ? 2, ? ( B1 ) ? 2 ? 1

所以G-S迭代法发散。
120

3.3 SOR迭代法(松弛法)
例题3-4 试判断以下简单迭代公式是否收敛:
1 (k ) 5 ? ( k ?1) x1 ? ? x2 ? ? ? 3 3 (1) ? 5 ? x ( k ?1) ? ? 1 x ( k ) ? 2 1 ? ? 2 2
( k ?1) (k ) (k ) ? x ? ? 0.7 x ? 0.5 x ? 1 1 2 ?4 (2) ? ( k ?1) (k ) (k ) x ? ? 0.1 x ? 0.2 x ? ? 2 1 2 ?3

解 (1)迭代矩阵

?1/ 3? ? 0 M ?? ?, ? 1/ 2 0 ? ?

M

?

? 1/ 2 ? 1

故迭代公式收敛。 (2)迭代矩阵
??0.7 0.5? M ?? , ? ? ?0.1 0.2?

M 1 ? max{0.8,0.7} ? 0.8 ? 1
121

故迭代公式收敛。

3.3 SOR迭代法(松弛法)
例题3-5 试判断简单迭代公式 是否收敛。
? 0 ?2? , 解:迭代矩阵 M ? ? ? ? ?3 0 ?
( k ?1) (k ) ? x ? ? 2 x ? 1 2 ?5 ? ( k ?1) (k ) x ? ? 3 x ?5 ? ? 2 1

? 2 ?I ? M ? ? ?2 ? 6 ? 0 3 ?

? ?? 6
故迭代发散。

? (M ) ? max{ 6, ? 6} ? 6 ? 1,

122

3.3 SOR迭代法(松弛法)
? 计算第k+1个分量xi(k+1)时将用GS方法计算的到的 结果记作 xi( k + 1) ,再引入参数ω,令:

x

( k + 1) i

= (1- w) x

(k ) i

+ wx

( k + 1) i

得到迭代公式

w ( k + 1) (k ) xi = (1- w) xi + (bi aii

( k + 1) a x 邋 ij j j= 1

i- 1

n j= i+ 1

aij x j (k ) )

(i = 1, 2,L , n)
123

3.3 SOR迭代法(松弛法)
矩阵表示 x(k+1)=(1-ω) x(k)+ω(Lx(k+1)+Ux(k) + g) 迭代矩阵: Bω=(I- ωL )-1((1-ω)I + ωU) ω=1即G-S法

124

3.3 SOR 收敛性
? 松弛法收敛的必要条件是0<ω<2
? 若系数矩阵具有严格对角优势,松弛因子ω满足 0<ω≤1,则松弛法收敛。 ? 若A是实对称矩阵,对角元为正,则当0<ω<2时 松弛法收敛的充分必要条件是A正定

125

对角占优阵
?严格对角占优矩阵
aii ? ? aij , i ? 1, 2,
j ?1 j ?i n

,n

126

3.3 SOR迭代法(松弛法)
简单迭代法收敛的充分条件:设方程组 Ax=b,构 造简单迭代公式为 x(k+1)= Mx(k)+g. 则: (1)若迭代矩阵M满足||M||1< 1 或 ||M||∞<1,则 简单迭代公式收敛。 (2)若系数矩阵A或其转臵AT是严格对角占优矩 阵,则雅可比迭代公式和高斯 -赛德尔迭代公式都 收敛。

(3)若系数矩阵A对称正定,则高斯-赛德尔迭代 公式收敛。
127

3.3 SOR迭代法(松弛法)
例3-5 试对以下方程组建立收敛的简单迭代公式:
?0.8 x1 ? 2 x2 ? 8 (1) ? ? x1 ? 0.5 x2 ? 4
? 4 x1 ? x2 ? 10 x3 ? 13 ? (2) ? 2 x1 ? 10 x2 ? x3 ? 11 ? 10 x ? 6 x ? 5 x ? 11 1 2 3 ?

解 (1) 将方程调整次序后得到等价的方程组
? x1 ? 0.5 x2 ? 4 ? ?0.8 x1 ? 2 x2 ? 8

其系数矩阵是严格对角占优的,构造雅可比迭代 公式一定收敛。迭代公式为
( k ?1) (k ) ? x ? 4 ? 0.5 x ? 1 2 ? ( k ?1) (k ) (k ) x ? (8 ? 0.8 x ) / 2 ? 4 ? 0.4 x ? ? 2 1 1

128

3.3 SOR迭代法(松弛法)
(2) 方程组
? 4 x1 ? x2 ? 10 x3 ? 13 ? ? 2 x1 ? 10 x2 ? x3 ? 11 ? 10 x ? 6 x ? 5 x ? 11 1 2 3 ?
?10 x1 ? 6 x2 ? 5 x3 ? 11 ? ? 2 x1 ? 10 x2 ? x3 ? 11 ? 4 x ? x ? 10 x ? 13 2 3 ? 1

调整方程的次序得等价的方程组 由此构造雅可比迭代公式

( k ?1) (k ) (k ) ? x1 ? (11 ? 6 x2 ? 5 x3 ) /10 ? ( k ?1) (k ) (k ) ? (11 ? 2 x1 ? x3 ) /10 ? x2 ? ( k ?1) (k ) (k ) x ? (13 ? 4 x ? x 1 2 ) /10 ? 3

迭代矩阵

?0.6 0.5 ? ? 0 ? ? M ? ? ?0.2 0 0.1 ? ? ?0.4 0.1 0 ? ? ?
129

由于||B||1=0.7<1,所以迭代公式是收敛的。

3.4 顺序高斯消去法
(消元过程按方程和未知量的自然顺序进行) 3.4.1 顺序高斯消去法举例 3.4.2 一般情况的计算过程

130

3.4 顺序高斯消去法
高斯消去法求解AX=b的基本思想: 先消元,即按一定的规律逐步消去未知量,将 方程组化为等价的上三角形方程组;然后进行回 代,即由上三角形方程组逐个求出

增广矩阵(A b)

初等行变换
131

3.4.1 顺序高斯消去法举例
例3-6 用顺序高斯消去 法求解方程组

解 消元过程:对增广矩阵[ A,b ]做消元计算,

132

3.4.1 顺序高斯消去法举例
[ A,b ]=
2 0 0 0 4 2 0 0 0 1 1

2 0 0 0
2 0.5 1.5 2 - 0.75 4.75 1 0.75 1.25

4 2 1 1
4 2
0 0

0 1 1 2 0.5 1.5 3 -0.5 5.5 2 1 2
0 2 1 0.5 1 1.5

2 0
0 0

2 -0.75 4.75 0 1.125 -1.125

133

3.4.1 顺序高斯消去法举例
消元过程结束,得上三角形方程组

134

3.4.1 顺序高斯消去法举例
回代过程:求解上 三角形方程组

所以方程组的解为
135

3.4.1 顺序高斯消去法的规则
? 基本的顺序高斯消去法,其消元过程遵守如下 规则:
?依从左到右、自上而下的次序将主对角元

下方的元素化为零。
?不作行交换,也不用非零数乘某行。 ?第 k列消元(将主对角元下方化零)时,将

下面各行分别减去第 k 行的适当倍数,不作 其他变换。
136

练习:
用顺序高斯消去 法求解方程组

2 1 5 3 1 6 0 1 2 2 3 0 -1/2 0 0 1 5 11 0 1 -1 1 9/2 -5 2 3 0 -1/2 0 -3/2 0 1 1 9/2

2 17/2

得上三角形方程 组,再回代得:
137

x1= –1 , x2= 1, x3= 5

3.4.2 一般情况的计算过程
常数项为 bk (k=1,2,…,n)

先看4阶的情形 设A=(aij)4非奇异,A(1)=A,b(1)=b,则原方程组为:

A(1)X= b(1)

138

3.4.2 一般情况的计算过程
第1 列消元(设a11(1)≠0),作

得到

139

3.4.2 一般情况的计算过程

其中



,则计算公式合写为

140

3.4.2 一般情况的计算过程
第2 列消元(设a22(2)≠0),作

化为

141

3.4.2 一般情况的计算过程

其中



,则计算公式合写为

142

3.4.2 一般情况的计算过程
第3 列消元(设a33(3)≠0),

化为上三角形方程组

其中

,计算公式写为

143

3.4.2 一般情况的计算过程
对于n 阶的一般情形

144

3.4.2 一般情况的计算过程
一般,对于n 阶的方程组,需要对前n - 1列进行 消元。设已进行了第1列至第k–1列消元,

145

3.4.2 一般情况的计算过程
第k列消元, 将第i 行 减去第k 行的l i k 倍,

(i=k+1,…,n)

146

3.4.2 一般情况的计算过程

147

3.4.2 一般情况的计算过程
第n-1列消元后

得到与原方程组同解的上三角形方程组

148

3.4.2 一般情况的计算过程
消元过程结束。进行回代过程,由上三角形方程 组依次求出xn ,xn-1 ,…,x1 。
若 则 一般由第k个方程得

149

顺序高斯消去法的计算步骤:
(1)消元过程。对 k=1,2,… ,n-1 计算 ①若akk=0,则算法失败, 结束计算;否则做下步。 ②对i=k+1,k+2,…,n计算

(2)回代过程。对 k=n,n-1,… ,1 计算

150

顺序高斯消去法的计算步骤:

称为消元过程的主元素。其中
151

3.4.2 顺序高斯消去法的适用条件
考察矩阵 的元 k阶顺序主子阵 Ak , 答: 问 题 1. A 消 过程能按 即A的左上角k行k 列的子阵,例如: 顺序进行到底的充要 1 2 3 A= 4 5 6 则A1=(1), A2= 1 2 , 条件是什么? A3=A

7 8 9

4 5

问题2.在原方程组的系 数矩阵中如何反映出这 个条件呢?

A的k阶顺序主子矩 阵Ak的行列式值是:

152

3.4.2 顺序高斯消去法的适用条件
定理 对于方程组(3.1),顺序消元过程能进行

到底的充要条件是系数矩阵A的顺序主子矩阵
Ak (k=1,2,…,n-1)非奇异。

153

顺序高斯消去法的数值稳定性如何?
例 用顺序高斯消去法解下面的 方程组,用3位有效数字计算 (假定机器字长为3)


误差: |ε(x1 )|=10,
-3-1000 2 |εr ( x1× )| =2=200% , =-2003≈-2000
误差太大!原因: 40-1000×40.1≈-40100 主元0.02绝对值比它 下面的元素20小得多。 结论:顺序高斯消 去法往往是数值不 稳定的。

回代求得 方程组的准确解是

由此产生选主元的思想。
154

3.5

选主元高斯消去法

基本思想: 每次消元前按一定的范围选取绝对值最大的元素 作为主元素。以便减少舍入误差的影响。 主要有列主元高斯消去法和全主元高斯消去法.

155

3.5.1 列主元高斯消去法
1、列主元高斯消去法:当k-1列消元后增广矩阵为

第k次消元时,先选列主元素,即在第k列 以下的各元素 中寻找绝对值最大的元素作为主元素 。即确定r,使 为

中之最大者。

156

3.5.1 列主元高斯消去法
只要det A≠0,就必能找到

若r>k 则交换第k行和第r行,然后进行消元。 例用列主元高斯消去法求解 方程组(用三位有效数字计算)
选主元



选主元
157

3.5.1 列主元高斯消去法

消元过程完成,得到上三角形方程组

再作回代可求得

158

3.5.1 列主元高斯消去法
练习: 用列主元高斯消 去法求解方程组 解
2 1 5 3 1 6 5 0 1 2 6 1 5 11 2 11 5 1 2 6 1 3 2 1 0 11 5 1

0 -1/5 3/5 14/5 0 3/5 -4/5 -17/5 5 6 2 11 0 3/5 -4/5 -17/5 0 0 1/3 5/3

5 6 2 11 0 3/5 -4/5 -17/5 0 -1/5 3/5 14/5

159

3.5.1 列主元高斯消去法
得上三 角形方程组

再回代得

x1= –1 , x2= 1, x3= 5
160

3.5.1 列主元高斯消去法
列主元高斯消去法(简称列主元消去法) 的计算步骤:

(1)消元过程。对 k=1,2,…,n-1 作下列运算: ①按列选主元: 确定 r,使满足

若ark=0,系数矩阵奇异,则停止计算并转结束;

②行交换: 若r>k , 则交换第r行和第k行,即对 j=k,k+1,…,n+1 计算
161

3.5.1 列主元高斯消去法
③消元计算:对 i=k+1,k+2, …,n, j=k+1,k+2, …,n+1 计算

(2)回代过程。若ann=0,则系数矩阵为奇异, 停止计算并转结束;否则计算

162

3.5.2 全主元高斯消去法
简称全主元消去法,在第k列(k=1,2, …,n-1)消元 时,从系数矩阵的右下角(n-k+1)阶子矩阵中,选 取绝对值最大的元素

不包含 常数列

163

3.5.2 全主元高斯消去法
例:用全主元高斯消 去法求解 方程组(用 三位有效数字计算)

选主元

选主元

164

3.5.2 全主元高斯消去法
得到等价的上三角形方程组

回代求得

165

3.5.2 全主元高斯消去法
对算法的说明

1. 系数矩阵为对称正定或是严格对角占优的情形 n阶矩阵A为对称正定的充要条件是AT=A且其一切 顺序主子式均大于0。
n阶矩阵A为严格对角占优矩阵是指其每个主对 角元的绝对值大于同一行其他元素绝对值之和, 即

一阶严格对角占优矩阵指一个非零数。
166

3.5.2 全主元高斯消去法
例如

由于第1行 第 2行
第 3行 所以B是严格对角占 优矩阵。

C不是严格对角占优矩阵
当方程组的系数矩 阵为对称正定或严格 对角占优矩阵时,可 用顺序高斯消去法求 解。
167

3.5.2 全主元高斯消去法
2.病态方程组 如果一个方程组的系数或者常数项只要有一点 微小的误差就会引起解的巨大变化,这样的方程 组就称为是病态的。


解就变成

由于右端系数有误差, 解的最大误差达200%。 方程组变成 对于病态方程组,前面的算法都可能失效, 168 需要使用一些改进精度的算法。

3.6.1 行列式的计算
⒈行列式的计算: p为消元过程中交换行(列)的次数。 例 列主元 法的消 元过程 计算过程有2次行交换,故p=2,主元为5,2.01, 1.78 det A ≈ (-1)2×5×2.01×1.78=18.0 (准确值det A=17.8325 )
169

3.6.1 行列式的计算
全主元 法的消 若A 元过程 计算过程进行了一次行交换,一次列交换,p=2,
则 det A≈(-1)2×5×2.10×1.70=17.9

170

3.6.2 逆矩阵的计算
当det A≠0时,A存在逆矩阵,设

求A-1就相当于求解n个线性方程组 s=1,2,…,n
171

3.6.2 逆矩阵的计算
具体计算时,令

同时求解n个方程组,并将解x(s) 存入原来 es 的位 臵(s=1,2,…,n)。 为了求逆方便,消元时可将主元上方也消为0,并 将第k行除以akk (k=1,2,…,n),使主元位臵化为1, 最后增广矩阵化为 这=A 种算法无须回 -1 代,称做高斯 - 约 当消去法。 172 P52-10

3.7 解三对角线性方程组的追赶法
三对角线性方程组

其系数矩阵为三对角矩阵
173

3.7 解三对角线性方程组的追赶法

当三对角矩阵满足条件

时,其所有顺序主子矩阵均非奇异。

174

3.7 解三对角线性方程组的追赶法
因此,三对角方程组的求解可以直接用顺序高 斯消元法而不必选主元。 将顺序消元过程略加修改:每次消元前将主元 素先化为1,即做“归一化”处理。

175

3.7 解三对角线性方程组的追赶法
第k次消元计算时增广矩阵第k、k+1两行形为

1 0
将第 k 行除以 b k , 第 k+1行减去第 k 行的 a k+1 倍。 计算公式:

176

3.7 解三对角线性方程组的追赶法
一般计算公式为 (令初值c n=0 ) (追)

消元结束 后,增广 矩阵化为
177

3.7 解三对角线性方程组的追赶法
消元结束后,增广矩阵化为

,得方程组

回代 计算 : 上述计算过程也称为追赶法。

(赶)

178

开始

输入 a, b, c, d d 1 ? d 1 / b1 ,c1 ? c1 / b1 对k ? 2 , 3 , ? , n计算

追赶法的计 算框图:
在计算机上一

般用四个一维数 组a[ ],b[ ],c[ ],d[ ] 分别存放三对角 线元素及方程组 右端常数向量, 其中令a[1]=0, c[n]=0。

? bk ? ?d k ?c ? k

? bk ? c k ? 1 a k ? ( d k ? d k ?1a k ) / bk ? c k / bk xn ? d n

x k ? d k ? c k x k ?1 k ? n ? 1, n ? 2 , ? ,1 输出 x
结束
179

练习
? 用追赶法求解三对角线性方程组

180

3.7 解三对角线性方程组的追赶法
追赶法的优点: (1)大大节省存储空间 只需用4n个单元。而一般的方程组需要n2+n个 单元。 (2)大大减少了计算量

一般的高斯消元过程需要O(n3) 次乘除法运算, 而“追”的过程只用4n次乘除运算.
181

3.8 三角分解法

?3.8.1 ?3.8.2

矩阵的三角分解 用三角分解法解方程组

184

3.8.1 矩阵的三角分解
· 初等行变换和初等矩阵
轾 aaa 犏 例如 A =21 犏 a2a a 3 2 2 犏 犏 3 2 31 臌 3a3a a
3 2 11 1 1

轾1 a1 2 1 3a 1 r3 - r2 5 犏 a2 2 2 3a 2 揪 揪 井 犏 1 犏 犏 a23 a -a 5 a a - a5 33 22 23 12 13臌
轾 1 0 0 犏 L 32 = 犏 0 1 0 犏 犏 0 - 5 1 臌

a a - 5

做初等矩阵 有

轾 轾 1 0 30 2 11 1a1a a 犏 犏 犏 L 32 A = 犏 0 1 30 a2a a 2 2 1 2 犏 犏 犏 犏 0 5 1 3 2 3 1 3a3a a 臌 臌

33

轾1 a1 2 1 3a 1 犏 = 犏1 a2 2 2 3a 2 犏 犏 a23 a -a 5 a a - a5 2 2 23 11 2 3臌

a a - 5
185

3.8.1 矩阵的三角分解
消元过程有
A
(k )

第i行 - 第k 行的lik 倍 ( k + 1) 揪 揪 揪 揪 井 A

lik = a

(k ) ik

a

(k ) kk



轾 1 犏 犏 O 犏 犏 犏 令 Lk = 犏 犏 犏 犏 犏 犏 犏 犏 臌

1 - lk + 1,k M M - lnk O L 1 O 1

? ?

第k 行 第 k + 1行

Lk A(k ) = A( k + 1)

186

3.8.1 矩阵的三角分解
一般消元过程依次将 a21,a31,…,an1,a32,a42,…,an2,…,an n-1化为0,相当于用 一系列初等矩阵L1,L2,…, Ln 依次左乘A,得到上三角 阵A(n) 。 ( n) L L L L L A= A ,
n- 1 3 2 1

A = (Ln- 1 L L3L2L1 ) A
令 L = (L n- 1 L L3L 2L1 )- 1 = L1 L 2 L L n- 2 L
- 1 - 1 - 1

- 1

( n)

- 1 n- 1

A = L1- 1L2- 1 L Ln- 2- 1Ln- 1- 1 A( n)
187

3.8.1 矩阵的三角分解
轾 1 犏 犏 - l21 1 L1 = 犏 犏M O 犏 犏 - ln1 L 臌
轾 1 犏 犏 l21 1 - 1 犏 L1 = 犏 M O 犏 犏 ln1 L 臌

1

1

轾 1 犏 犏 l21 1 - 1 - 1 - 1 - 1 L1 L 2 L L n- 2 L n- 1 = 犏 犏 M O 犏 犏 ln1 ln 2 L 臌

1
188

3.8.1 矩阵的三角分解
得 A = LA( n)
轾 1 犏 犏 l21 1 犏 = 犏 M O 犏 犏 ln1 ln 2 L 臌
(1) 轾 a 11 犏 犏 犏 犏 犏 1 犏 犏 臌

a12 (1) L
(2) a22

a1n (1)
(2) a2 n M (n) ann

L O

定义 设A为n阶矩阵,若A有分解式 A=LU 其中L为下三角阵,U为上三角阵,则称之为A 的LU分解或三角分解。

三角分解
189

3.8.1 矩阵的三角分解
u11? u12 ? 0 u? ??1 ? 1 0 ?? 2 ? ?3 1n ??? 2 ? 3 ? ? ?1 ?2 3 ?? ?? ? ? ? ? 例如A= ? , ? ? ? ? ?? ? l 1 u ? u ? ? ? ? ? ? 5 ?4 1 ? 杜利特尔 ? ? 2n ? 0 ? 3? ? 5 ? 22 ? ? 2 3? ? ? ??21 ? 2 1 ?? 0 ?? ? ?? ? ?? ? ? ? 分解 ? ?? ? l l ? 1 u 这里有A的两种不同的三角分解,类似可举出很 nn ? (Doolittle) ? n1 n 2 ??

常用的两种三角分解:

多,一般,若A=LU是一个三角分解,任取与A ? l11 ? ?1 u12 ? u1n ? 克劳特分 同阶的非奇异对角矩阵 D ,则 ?l ? ? ? l 22 1 ? u2 n ? 21 ? ? ? ? -1U)=(LD)(D-1U)=L 解 A= (LDD ?? ?? ? ? ? ? 1U1 ? ?? ? (Crout) 也是A的三角分解。 l l ? l 1 ? nn ? ? ? n1 n 2
问:任一矩阵是否存在三角分解?若存在,是否惟一? 答:不一定存在。若存在,则不惟一。
190

3.8.1 矩阵的三角分解
定理3-6 n阶矩阵A存在唯一的杜利特尔分解和克 劳特分解的充要条件是A的顺序主子矩阵Ak (k=1,2,…,n-1)非奇异。 P52-12 推论 若定理3-6的条件满足,则A存在唯一的 LDU分解式,其中为L单位下三角阵,U为单位上 三角阵,D为对角阵。 即 A=LDU
轾 1 犏 犏 l21 1 犏 = 犏 M O 犏 犏 ln1 ln 2 L 臌 轾 d1 犏 犏 d2 犏 犏 O 犏 1 犏 臌 轾 1 u12 L 犏 犏 1 犏 犏 O 犏 dn 犏 臌 u1n u2 n M 1
191

反之也真。

3.8.1 矩阵的三角分解
?1 ?l 21 ?? ? ? ? ? l n1 1 ln2 ? ? u11 ?? ?? ? ?? ?? ? 1? ? u12 u22 ? u1n ? 杜利特尔 ? u2 n ? ? ? ? ? 分解 ? unn ? (Doolittle)

? l11 ?l 21 ?? ?? ? ? l n1

l 22 ln2

? ?1 u12 ? u1n ? ?? ? 1 ? u 2n ? ?? ?? ? ? ? ? ?? ? ? l nn ? ? 1 ?

克劳特分 解 (Crout)
u1n u2 n M 1
192

A=LDU

轾 1 犏 犏 l21 1 犏 = 犏 M O 犏 犏 ln1 ln 2 L 臌

轾 d1 犏 犏 d2 犏 犏 O 犏 1 犏 臌

轾 1 u12 L 犏 犏 1 犏 犏 O 犏 dn 犏 臌

3.8.1 矩阵的三角分解
例 试求矩阵
轾 2 4 A= 犏 犏 - 4 - 5 臌

的杜利特尔分解,克劳特分解,LDU分解。 解 设A的杜利特尔分解为
轾 轾 1 0 轾 u11 u12 2 4 犏 犏 = 犏 犏 l21 1 犏 0 u22 - 4 - 5 犏 臌 臌 臌

右端矩阵相乘,两端比较可依次求得 u11 ? 2, u12 ? 4 l 21 ? ?4 / u11 ? ?2 l21u11 = - 4, l21u12 + u22 = - 5,
u22 ? ?5 ? l21u12 ? ?5 ? (?2)4 ? 3
193

3.8.1 矩阵的三角分解
得杜利特尔分解为
轾 1 0轾 2 4 犏 A= 犏 犏 - 2 1犏 0 3 臌 臌

轾 2 4 由此可得LDU分解为 A= 犏 犏 - 4 -轾 5 臌 1
克劳特分解为

0 A= 犏 犏 - 2 1 臌

轾 2 0轾 1 2 犏 犏 犏 0 3犏 0 1 臌 臌

轾 2 0轾 1 2 犏 A= 犏 犏 - 4 3犏 0 1 臌 臌

P52-13

194

3.8.2 用三角分解法解方程组
1. 思想: 设A非奇异,并有三角分解A=LU,则方程组 Ax=b 就化为 LUx=b 只须求解两个简单的三角形方程组: (1)解Ly=b 求出 y ,(2) 解Ux=y ,求出x. 2.例 用杜利特尔分解法求解 ì ? ? ? ? 方程组 í ? 解
2 x1 + x2 - x3 = - 1 4 x1 - x2 + 3 x3 = 7

? ? ? ? 6 x1 + 9 x2 - x3 = - 3

首先求分解式。设系数矩阵的杜利特尔分解 为A=LU,即
195

3.8.2 用三角分解法解方程组
?2 1 ? 1? ? 1 ?4 ? 1 3 ? ? ? l ? ? ? 21 ? ?6 9 ? 1? ? ? ?l 31 0 1 l 32 0? ?u11 ?0 0? ?? 1? ?? ?0 u12 u22 0 u13 ? u23 ? ? u33 ? ?

右端两矩阵相乘, 并比较等式两端。

u13 ? ?1 u12 ? 1 u11 ? 2 l21 ?u11 4, ? l21 ? 4 / u11 ? 2 l31 ?u11 6, ? l31 ? 6 / u11 ? 3
- 1 = l21u12 + u22 ,
3 = l21u13 + u23 ,
l31u12 + l32u22 = 9,

u22 ? ?1 ? l 21 u12 ? ?1 ? 2 ? 1 ? ?3
u23 ? 3 ? l21u13 ? 3 ? 2 ? (?1) ? 5

l32 ? (9 ? l31u12 ) / u22 ? ?2
u33 ? ?1 ? l31u13 ? l32 u23 ? 12
196

- 1 = l31u13 + l32u23 + u33 ,

3.8.2 用三角分解法解方程组
轾 1 0 0 犏 L= 犏 2 1 0 犏 犏 3 - 2 1 臌 , 轾 2 1 - 1 犏 U= 犏 0 - 3 5 犏 犏 0 0 12 臌

其次,求解两个三角形方程组,
ì y1 =- 1 ? ? ? ? = 7 í 2 y1 + y2 ? ? ? ? ? 3 y1 2 y2 + y3 = - 3 ì 2 x1 + x2 - x3 = y1 ? ? ? ? í - 3 x2 + 5 x3 = y2 ? ? ? 12 x3 = y3 ? ?

Ly ? b, Ux ? y



y1 = - 1 y2 = 9 y3 = 18 x1 = 1/ 2 x2 = - 1/ 2 x3 = 3/ 2
197

3.8.2 用三角分解法解方程组
3 一般计算公式:设A的杜利特尔分解式为
轾 a11 犏 犏 a21 犏 犏 犏 犏 an1 臌 a12 L L L ann a1n a2 n a22 L an 2 L 轾 1 犏 犏 l21 1 犏 = 犏 M M O 犏 犏 ln1 ln 2 L 臌 轾 u11 u12 L 犏 犏 u22 L 犏 犏 O 犏 1 犏 臌 u1n u2 n M unn

比较两端第一行,得U的第一行元素
u1 j ? a1 j , j ? 1,2,?, n

比较两端第一列有ai1 = li1u11 ,得L的第一列元素

li 1 ? ai 1 / u11 ,

i ? 2,3,?, n
198

若U的前k-1行、L的前k-1列已求,

3.8.2 用三角分解法解方程组
欲求U的第k行元素 ukj ( j ? k ) ,考虑
骣 a11 a12 珑 珑 珑 a21 a22 珑 珑 珑 珑 L L 珑 珑 珑 ak1 ak 2 珑 珑 珑 珑 L L 珑 珑 珑 珑 an1 an 2 桫 L L L L L L a1 j L a2 j L L akj L anj L L L L a1n 鼢 鼢 鼢 a2 n 鼢 鼢 鼢 鼢 L 鼢 鼢 鼢 = 鼢 akn 鼢 鼢 鼢 ÷ L ÷ ÷ ÷ ÷ ann ÷ 骣 1 l21 L 0 1 L L L L 0 0 L L L L 骣 u11 u12 ? ? ? 0 u22 ? ? ? L ? L L ? ? 0 ? 0 0 ? ? ÷? L ÷? L L ÷ ? ÷ ? ÷ ? 1 ÷? 0 0 桫 0 0 L L L L L L u1 j L u2 j L L ukj L 0 L L L L

lk 1 lk 2 L L L L ln1 ln 2 L 桫

lkk = 1 L L L lnj L

u1n ÷ ÷ ÷ u2 n ÷ ÷ ÷ ÷ L ÷ ÷ ÷ ÷ ukn ÷ ÷ ÷ ÷ L ÷ ÷ ÷ ÷ unn ÷

akj = lk1u1 j + lk 2u2 j + L + lk k- 1uk- 1 j + ukj
于是
uk j ? ak j ? ? l k s us j , j ? k , k ? 1, ? , n
s ?1
199

k ?1

3.8.2 用三角分解法解方程组
骣 a11 珑 珑 珑 a21 珑 珑 珑 珑 L 珑 珑 珑 ai1 珑 珑 珑 珑 L 珑 珑 珑 珑 an1 桫 a12 L a22 L L ai 2 L L L L

欲求L的k列元素 lik (i > k ) ,考虑
a1k a2 k L aik L ank L L L L L L a1n 鼢 鼢 鼢 a2 n 鼢 鼢 鼢 鼢 L 鼢 鼢 鼢 = 鼢 ain 鼢 鼢 鼢 鼢 L 鼢 鼢 ÷ ÷ ann ÷ 骣 1 l21 L li1 L 0 1 L L 0 0 L L 0 0 L li 2 L L L L L lik L lnk L L L L

an 2 L

ln1 ln 2 L 桫

骣 u11 u12 ? ? ? 0 u22 ? ? ? L ? L L ? ? 0 ? 0 0 ? ? ? L ? L L ? 鼢 ? 鼢 ? ? 1鼢 0 0 桫

L L L L L L

u1k u2 k L ukk L 0

L L L L L L

u1n ÷ ÷ ÷ u2 n ÷ ÷ ÷ ÷ L ÷ ÷ ÷ ÷ ukn ÷ ÷ ÷ ÷ L ÷ ÷ unn

aik = li1u1 j + li 2u2k L + li k- 1uk- 1k + lik ukk
于是 l i k ? (ai k ? ? l i s us k ) / uk k , i ? k ? 1, k ? 2, ?, n
s ?1
200

k ?1

3.8.2 用三角分解法解方程组
求得L、U后再求解方程组Ly=b, Ux=y
ì y1 = b1 y1 = b1 , y2 = b2 - l21 y1 ,L ? ? ? ? l21 y1 + y2 = b2 ? ? ? k ?1 ............... ? ? y ? b ? l y , í k k k j j lk1 y1 + lk 2 y2 + L + lk k - 1 yk - 1 + yk = bk ? ? j ?1 ? ? ............... ? ? k ? 1,2,? , n ? ? ? ? ln1 y1 + ln 2 y2 + L + ln n- 1 yn- 1 + yn = bn

?

201

3.8.2 用三角分解法解方程组
ì u11 x1 + u12 x2 + L + u1n xn = y1 ? ? ? ? u22 x2 + L + u2 n xn = y2 ? í ? LLLL ? ? ? unn xn = yn ? ? ?

x k ? ( yk ?

j ?k ?1

?u

n

kj

x j ) / uk k
202

3.8.2 用三角分解法解方程组
用杜利特尔分解法解方程组的计算步骤: (1)三角分解。对k=1,2,…, n,计算 u ? ?1 ? ?u
11 12

u1n ? ? ?l ? ? k ?1 ? 设A=LU= u ? u 1 22 2n ? 21 ? ? ? u ? a ? l u , j ? k , k ? 1 , ? , n k j k j k s s j ? ?? ? ? ? ? ? ? ? ? s ?1 ? ? ? ? ? k ?1 unn ? l n 1 l n 2 ? 1? ? ?l ? ( a ? ? l u )/ u , i ? k ? 1, k ? 2, ? , n

?
s ?1

? ?

ik

ik

?

is

sk

kk

(2)求解。 ①解Ly=b。对k=1,2,…,n 计算 ②解Ux=y。对k=n,n-1,…,1 计算
x k ? ( yk ?
j ?k ?1

yk ? bk ? ? l k j y j
j ?1

k ?1

?u

n

kj

x j ) / uk k

203

开 始

输入A , b
k ?1

对 j ? k , k ? 1,? , n ak
j

计算

? ak

j

?

k ?1 s ?1

? ak sas j
k ?1 s ?1

对 i ? k ? 1, k ? 2,? , n 计算

ai k ? ( ai k ? ? ai s a s k ) a k k

bk ? bk ?

k ?1 j ?1

? ak jb j
T
k ? k ?1

k ?n F
n

紧 凑 存 贮 的 杜 利 特 尔 分 解 法

对 k ? n, n ? 1,? ,1 计算 bk ? (bk ?
j ? k ?1

? ak jb j ) ak k

输出 b

204

3.9 线性方程组的最小二乘法
? 线性方程组AX=b 无解,称之为不相容方程组或 矛盾方程组。 ? 方程个数> 未知数个数,称为超定方程组。

n> m

205

3.9 线性方程组的最小二乘法
n> m

如果有向量X

则称向量X为AX=b的最小二乘解
206

3.9 线性方程组的最小二乘法
定理3-7 线性方程组AX=b必存在最小二乘解,且 即为方程组ATAX=ATb 的解,当A列满秩时, AAT是非奇异的,这时有唯一解。 证明:设方程组ATAX=ATb 的一个解为C,则有 ATAC=ATb ,设u是任一个m维的向量,且u=c+z 则:

207

3.9 线性方程组的最小二乘法
? 求 线性方程组的最小二乘法的步骤:

1、确定系数矩阵A和AT。计算ATA和ATb
2、求方程组ATAX=ATb的解

方程组ATAX=ATb 叫做超定方程组 AAX=b的正则方程 组或法方程组
208

3.10 病态方程组及条件数
如果一个方程组的系数或者常数项只要有一点微 小的误差就会引起解的巨大变化,这样的方程组 就称为是病态的。 例
解就变成

由于端常数项有误差, 解的最大误差达200%。 方程组变成 对于病态方程组,前面的算法都可能失效, 需要使用一些改进精度的算法。 209

3.10 病态方程组及条件数
? 条件数 ? Cond(A)=║A-1║║A║ ? Cond(A)≥1 ? Cond(cA)=Cond(A),c≠0 ? 扰动分析表明:条件数不大, 扰动对解的影响不大;条 件数越大,扰动对解的影响也越大.这就是说条件数是 方程组敏感性以及病态或良态的度量. ? 上例中Cond (A)1=2.00012×104

210

211

第五章 插值法和曲线拟合
主讲:孙剑

聊城大学计算机学院

本章主要内容:
?

5.1 5.3 5.4 5.5 5.6

插值法的基本理论 牛顿均差插值多项式 三次Hermite插值 三次样条插值 B样条插值
213

?
? ? ? ?

5.2 拉格朗日插值多项式

?

5.7 曲线拟合的最小二乘法

5.1 插值法的基本理论
?

5.1.1 插值问题及代数多项式插值

?

5.1.2 插值多项式的误差

214

5.1.1 插值问题及代数多项式插值
?函数y=f(x)给出一组函数值 yi = f ( xi ), i = 0,1,L , n
x: y: x0 x1 x2 …… xn y0 y1 y2 …… yn
插值基点

插值区间

其中,x0 , x1 , x2 ,…, xn是区间[a,b]上的互异点,要 构造一个简单的函数 j ( x ) 作为f(x)的近似表达式, 使满足 j ( xi ) = yi , i = 0,1,L , n 若x∈[a,b],需要计算f(x) 的近似值 j ( x ) ,则称x 插值条件 被插值函数 插值函数 为插值点。 插值原则
215

5.1.1 插值问题及代数多项式插值
y0

y1



yn

y = j (x )

y=f(x)

x0 x1 …

xn

图5-1 插值问题
216

5.1.1 插值问题及代数多项式插值
?当选择代数多项式作为插值函数时,称为代数多 项式插值问题。 设函数y=f(x)在[a,b]有定义, 且已知在n+1个点 a≤x0<x1<……<xn≤b上的函数值y0, y1,……,yn.,要求 一个次数不高于n的多项式

Pn ( x) ? a0 ? a1 x ? a2 x ? ?? an x
2

n

使满足插值原则 Pn ( xi ) ? yi , i ? 0,1,?, n 称Pn(x)为f(x)的n次插值多项式。

??这样的插值多项式是否存在、唯一? ?

217

5.1.1 插值问题及代数多项式插值
定理5.1 在n+1个互异基点处满足插值原则且次数 不超过n的多项式Pn(x)是存在并唯一的。 证明: 由 Pn ( xi ) = yi

,

i = 0,1,L , n



2 n ì ? a + a x + a x + L + a x 0 1 0 2 0 n 0 = y0 ? ? 2 n ? a + a x + a x + L + a x ? 0 1 1 2 1 n 1 = y1 ? í ? LL ? ? ? 2 n ? a + a x + a x + L + a x 1 n 2 n n n = yn ? ? 0
218

5.1.1 插值问题及代数多项式插值
证明:
其系数行列式 1 x x 2 K x n 0 0 0
V ( x0 , x1 ,L , xn ) =

范德蒙 行列式
=

1 x1 1 xn

x K x LL
2 n

2 1

n 1

0? j i n

?

( xi - x j )

0

x K x

n n

因此方程组存在唯一的解 a0 , a1 ,L , an 因此Pn(x)存在并唯一。
219

5.1.1 插值问题及代数多项式插值
定理5.1 在n+1个互异基点处满足插值原则 且次数不超过 n的多项式 Pn(x)是存在并唯一

的。

220

5.1.2 插值多项式的误差
Rn ( x) ? f ( x) ? Pn ( x) 也称为P (x)的余项。 截断误差: n

定理5.2 设函数f(x)在包含基点x0 , x1 ,…, xn的区间 [a,b]上连续,在(a,b)上具有 n+1阶导数,Pn(x)为满 足插值原则的n次插值多项式,则对任一点x∈[a,b], 总存在相应的点ξ∈(a,b),使 f ( n?1) ?? ? Rn ( x ) ? ? n?1 ( x )
( n ? 1)!

其中 w ( x) = n+ 1

? (xi= 0

n

xi ) = ( x - x0 )( x - x1 ) L ( x - xn )

221

5.1.2 插值多项式的误差
证明: ①当x=xi (i=0,1,…,n) 时,由插值条件知Rn(xi)=0, 故结论显然成立。

( n ?1) ? f ?? ②当x是[a,b]上任一个固定点,但不是插值基点时, Rn ( x ) ? ? n?1 ( x ) ( n ? 1)!

作辅助函数

Rn ( x) (t ) = ? f (t( )x -- P () t) (1t) )L ( x - xn ) x =- ( x - x0 )( w xx 其中 wn+F ni n+ 1 1(x wn+ 1 ( x) i= 0
n

F(t)在区间[a,b]上至少有n+2个互异的零点 x, x0 , x1 ,…, xn

222

5.1.2 插值多项式的误差

根据罗尔定理,在F(t)的两个零点之间至少有一点 使F(t)的一阶导数为0,即 F?(t)在(a,b)内至少有n+1个互异的零点, F??(t) 在(a,b)内至少有n个互异零点, 依此类推,可知F(n+1)(t)在(a,b)内至少有一个零点ξ, 则 n R ( x ) n R ( x ) (n + 1) ( n+P 1) (t ) n F t ) = f ( t ) (= tx )n0 w ( t ) = ( t x ) = ( t x )w L (t1) Fn( ( x ) = f ( x ) (1n +n+ 1)! ? +1 in 0 )(t - x w ()x) w i= 0 n1 +( 1x n+
f ( n?1) ?? ? Rn ( x ) ? ? n?1 ( x ) ( n ? 1)!

从而结论成立。
223

5.1.2 插值多项式的误差

定理5.2 设函数f(x)在包含基点x0 , x1 ,…, xn的区间 [a,b]上连续,在(a,b)上具有n+1阶导数,Pn(x)为满 足插值原则的n次插值多项式,则对任一点x∈[a,b], 总存在相应的点ξ∈(a,b),使

? f ?? Rn ( x ) ? ? n?1 ( x ) ( n ? 1)!
( n ?1)

其中 wn+ 1 ( x) =

? (xi= 0

n

xi ) = ( x - x0 )( x - x1 ) L ( x - xn )
224

5.1.2 插值多项式的误差
推论 当f(x)是次数不超过n的多项式时,其n次插 值多项式就是f(x)本身。 证明: 因为f(n+1)(x)=0, f ( n?1) ?? ?
( n ? 1)!

Rn ( x ) ?

? n?1 ( x )

从而Rn(x)≡0,

Pn(x)≡f(x)。
225

5.1.2 插值多项式的误差
f ( n?1) ?? ? 误差公式的用法: Rn ( x ) ? ? n?1 ( x ) ( n ? 1)!

如果

max f ( n+ 1) (t ) ? M ,则截断误差估计为
a# t b

M ? n ?1 ( x ) Rn ( x ) ? ( n ? 1)!

其中, M为 f ( n+ 1) (t ) 在区间(a,b)内的极大值 和两个端点处的值的最大值。
226

5.1.2 插值多项式的误差
设min{ x0 , x1 ,…, xn}=a,max{x0 , x1 ,…, xn}=b

当插值点x∈(a,b)时称为内插,否则称为外插。
y1 … yn



y=ψ(x)
y=f(x)

y0

x0 x1 x

xn x
227

5.1.2 插值多项式的误差
例5-1 若利用 x 在100和144的值,运用多项式插 值计算 110 的值,估计误差是多少? 解: 插值基点x0=100, x1=144, 插值点 110.
f ( 2 ) ?? ? f ( 2 ) ?? ? R1 ( x ) ? ? 2 ( x) ? ( x ? 100)( x ? 144) 2! 2! 3 1 1 1 -1 2 2 f ( x) = x , f '( x) = x , f ''( x) = - x = 2 4x x 1 14 (2)

100# x 144

max

f

(x ) =

100# x 144

max

4x x

=

4000

1 1 | ( x 100)( x - 144) | | R1 ( x) |? ? 2 4000 1 | (110 - 100)(110 - 144) | = 17 = 0.0425 | R1 (110) |? 8000 400

228

5.2 拉格朗日插值多项式

?

5.2.1 线性插值

?
?

5.2.2 二次插值
5.2.3 n次拉格朗日插值

229

5.2.1 线性插值
1 线性插值----n=1时的代数多项式插值 x x0 x1 已知f(x0)=y0,f(x1)=y1, x0≠x1 y y0 y1 要构造线性函数 P1(x)=a0 + a1 x , 使满足插值条件 P1(x0)=y0 , P1(x1)=y1 .
y -? y 11 0 y y 0 (x - x ) P ( x ) y = P1 ( x ) ? 0 y0 ? (x ? x 1 0 0 ) (线性插值多项式) x -? x x x 11 00

x - x0 x - x1 P + y1 1 ( x) = y0 x0 - x1 x1 - x0

拉格朗日线性插值多项式 230

5.2.1 线性插值
x ? x0 x ? x1 P1 ( x ) ? y0 ? y1 ? y0 l0 ( x ) ? y1l1 ( x ) x0 ? x1 x1 ? x0

拉格朗日线性插值多项式 公式的结构:它是两个一次函数的线性组合:

x ? x1 l0 ( x ) ? , x0 ? x1
基函数的性质

x ? x0 l1 ( x ) ? x1 ? x0
,

线性插值 基函数

li ( xi ) = 1 , li ( x j ) = 0 ( j ? i)

i, j 0,1

231

5.2.1 线性插值
线性插值的几何意义:用通过点 A( x0 , f ( x0 )) 和 B( x1, f ( x1 )) 的直线近似地代替曲线y=f(x)
y=f(x) p(x)=ax+b A(x.0,f(x.0)) B(x.1,f(x.1))

由解析几何知道,这条直线用点斜式表示为
y1 - y0 p( x) = y0 + ( x - x0 ) x1 - x0
x - x0 x - x1 p( x) = y0 + y1 x0 - x1 x1 - x0
232

5.2.1 线性插值
2、线性插值多项式的截断误差为 :
f ' ' (? ) R1 ( x ) ? f ( x ) ? P1 ( x ) ? ( x ? x0 )( x ? x1 ) 2!

ξ是在包含x, x0 , x1 的区间内某数。

233

5.2.1 线性插值
例5-2 给定函数y=lnx在两点10、11的值如下表, 用线性插值求ln10.5的近似值,并估计截断误差。 x 10 11 解: f(x)=ln x ,x0=10, x1=11, x=10.5 y 2.303 2.398
x ? x0 x ? x1 P1 ( x ) ? y0 ? y1 x0 ? x1 x1 ? x0
10.5 - 10 10.5 - 11 + 2.398 = 2.350 5 2.303? 11- 10 10 - 11
10# x 11

ln 10.5≈P1(10.5)=
f ? ( x) = 1 , x

1 f "( x) = - 2 x

max f "(x ) ? 1/102

0.01
234

0.01 R1 (10.5) ? (10.5 10)(10.5 - 11) = 0.001 25 2!

5.2.1 线性插值
x ? x0 x ? x1 P1 ( x ) ? y0 ? y1 ? y0 l0 ( x ) ? y1l1 ( x ) x0 ? x1 x1 ? x0

(回顾)

拉格朗日线性插值多项式 公式的结构:它是两个一次函数的线性组合:

x ? x1 l0 ( x ) ? , x0 ? x1
基函数的性质

x ? x0 l1 ( x ) ? x1 ? x0
,

线性插值 基函数

li ( xi ) = 1 , li ( x j ) = 0 ( j ? i)

i, j 0,1

235

5.2.2 二次插值
二次插值----n=2时的代数多项式插值 已知f(x)在三个互异点x0 ,x1 ,x2的函数值y0 ,y1 ,y2
x y x0 y0 x1 y1 x2 y2
2

要构造次数不超过二次的多项式

P2 ( x) = a0 + a1x + a2 x
使满足插值条件

P2 ( xi ) = yi , (i = 0,1, 2)
236

5.2.2 二次插值
?公式的构造: 先对每个基点xi构造二次插值基函 数li(i=0,1,2),使满足
li (xi ) ? 1 , li (x j ) ? 0(j ? i) , i、j ? 0,1 ,2
( x ? x1 )( x ? x2 ) l0 ( x ) ? ( x0 ? x1 )( x0 ? x2 )

( x ? x0 )( x ? x2 ) l1 ( x ) ? ( x1 ? x0 )( x1 ? x2 )
( x ? x0 )( x ? x1 ) l2 ( x) ? ( x2 ? x0 )( x2 ? x1 )

237

5.2.2 二次插值
P2 ( x) ? y0l0 ( x) ? y1l1 ( x) ? y2l2 ( x)
? y0 ( x ? x0 )( x ? x2 ) ( x ? x0 )( x ? x1 ) ( x ? x1 )( x ? x2 ) ? y1 ? y2 ( x0 ? x1 )( x0 ? x2 ) ( x1 ? x0 )( x1 ? x2 ) ( x2 ? x0 )( x2 ? x1 )

拉格朗日二次插值多项式 P2(x)的截断误差为
f ' ' ' (? ) R2 ( x ) ? f ( x ) ? P2 ( x ) ? ( x ? x0 )( x ? x1 )( x ? x2 ) 3!

ξ是包含x0 , x1 , x2 , x的区间内某数。

238

5.2.2 二次插值
例5-3 已知函数y=f(x)的观测数据如表所示,试求 其拉格朗日插值多项式,并计算f(1.5)的近似值。
x 解: y 0 1 2 -1 2 4

( x - 1)( x - 2) ( x - 0)( x - 1) ( x - 0)( x - 2) P2 (P x) ( = 2 + (y - 1) + 4 y l ( x ) ? l ( x ) ? y l ( x ) 2 x) ? 0 0 1 1 2 2 (0 - 1)(0 - 2) (1- 0)(1- 2) (2 - 0)(2 - 1)
2x )( x ? x ) ( x ? 4 x -1 7 x +2 2 ? y ( x ? x0 )( x ? x2 ) ? y ( x ? x0 )( x ? x1 ) ? y= 0 1 2 ( x0 ? x1 )( x0 ? x2 ) ( x1 ? x0 )( x1 ? x2 ) ( x2 ? x0 )( x2 ? x1 )

f (1.5) ? P 2 (1.5)

4? 1.5

2

7? 1.5 2 = 0.5

239

5.2.2 二次插值

二次插值的几何意义:用经过3个点( x 0 , y0 ),( x 1 , y1 ),( x 2 , y2 ) 的抛物线 y = P( x) 近似代替曲线 y = f ( x) ,如下图所 示。因此也称之为抛物插值。

240

5.2.3 n次插值
作基点xi的n次插值基函数 ( i=0,1,…,n):
( x ? x0 )( x ? x1 )?( x ? xi ?1 )( x ? xi ?1 )?( x ? xn ) li ( x ) ? ( xi ? x0 )( xi ? x1 )?( xi ? xi ?1 )( xi ? xi ?1 )?( xi ? xn ) ??
n

x ? xj xi ? x j

j ?0 j?i

i ? 0,1,? , n

基函数具有性质: li (xi ) ? 1 , li (x j ) ? 0(j ? i) , i、j ? 0,1,?, n n次拉格朗日插值多项式
Pn ( x) ? y0l0 ( x) ? y1l1 ( x) ? ? ? ynln ( x) ? ? yi ?
i ?0 j ?0 j?i n n

x ? xj xi ? x j
241

5.2.3 n次插值

242

拉 格 朗 日 插 值 算 法 实 现

243

5.2.3 n次拉格朗日插值
作基点xi的n次插值基函数 ( i=0,1,…,n):
( x ? x0 )( x ? x1 )?( x ? xi ?1 )( x ? xi ?1 )?( x ? xn ) li ( x ) ? ( xi ? x0 )( xi ? x1 )?( xi ? xi ?1 )( xi ? xi ?1 )?( xi ? xn ) ??
n

x ? xj xi ? x j

j ?0 j?i

i ? 0,1,? , n

基函数具有性质: li (xi ) = 1 , li (x j ) = 0 n次拉格朗日插值多项式

(j ? i),
n

i、j
n

0,1,L , n
x ? xj xi ? x j
244

Pn ( x) ? y0l0 ( x) ? y1l1 ( x) ? ? ? ynln ( x) ? ? yi ?
i ?0 j ?0 j?i

5.3 牛顿均差插值多项式
?
? ?

5.3.1 均差概念及计算
5.3.2 牛顿型插值多项式 5.3.3 差分及等距基点的牛顿插值公式

245

5.3.1 均差概念及计算

? 拉格朗日插值多项式结构对称,使用方便。但由 于是用基函数构成的插值,这样要增加一个节点 时,所有的基函数必须全部重新计算,不具备递 推性,还造成计算量的浪费。这就启发我们去构 造一种具有递推性的插值多项式来克服这个缺点, 也就是说,每增加一个节点时,只需增加相应的 一项即可。这就是牛顿插值多项式。

246

5.3.1 均差概念及计算
目的:构造具有如下形式的插值多项式:

P n ( x) = a0 + a1 ( x - x0 ) + a2 ( x - x0 )( x - x1 ) + L + an ( x - x0 )( x - x1 ) L ( x - xn- 1 )
它满足递推性:
P n ( x) = P n- 1 ( x) + an ( x - x0 )( x - x1 ) L ( x - xn- 1 )

247

5.3.1 均差概念及计算
一阶均差

x j ? xi (一阶差商 ) f [ x j , xk ] ? f [ xi , x j ] 二阶均差 f [ xi , x j , xk ] ? xk ? xi
例如设
f (0) = 1, f (2) = 5, f (4) = 9

f [ xi , x j ] ?

f ( x j ) ? f ( xi )

则 f [0, 2] = ( f (2) - f (0)) (2 - 0) = 2

f [2, 4] = ( f (4) - f (2)) (4 - 2) = 2 f [0, 2, 4] = ( f [2, 4] - f [0, 2]) (4 - 0) = 0

248

5.3.1 均差概念及计算
k阶均差
f [ xi ?k ?1 , ?, xi ?1 , xi ] ? f [ xi ?k , xi ?k ?1 , ?, xi ?1 ] f [ x i ? k , x i ? k ?1 , ? , x i ] ? xi ? xi ?k

例如

零阶均差定义为函数值本身,即 f [ xi ] ? f ( xi ) 均差具有对称性:任意改变基点的次序后其值不 变。例如 f [0,2,4] = f [2,0,4] = f [4,2,0]等等。
249

f [2, 4,5] - f [0, 2, 4] f [0, 2, 4,5] = =- 1 5- 0

5.3.1 均差概念及计算 粗线框出的部分在计算机
均差表1
xi f(xk) 1阶 2阶

上可存入二维数组
3阶
f [ x 2 , x 3 , x4 ] ?

4阶
f [ x 3 , x4 ] ? f [ x 2 , x 3 ] x4 ? x 2

x0 f(x0)
x1 f(x1) f [x0,x1 ] x2 f(x2) f [x1,x2] f [x0,x1,x2]

f [ x2 , x3 , x4 ] ? f [ x1 , x2 , x3 ] f [ x1 , x2 , x3 , x4 ] ? x4 ? x1

x3 f(x3) f [x2,x3 ] f [x1,x2,x3 ] f [x0,x1,x2,x3 ]
x4 f(x4) f [x3,x4] f [x2,x3,x4 ] f [x1,x2,x3,x4 ] f [x0,x1,x2,x3,x4 ]

┊ ┊








250

5.3.1 均差概念及计算
均差表2
xk f ( xk ) 一阶差商 x0 f ( x0 )
f [ x0 , x1 ]

二阶差商

三阶差商

四阶差商

x1 f ( x1 )
f [ x1 , x2 ]

f [ x0 , x1 , x2 ] f [ x0 , x1 , x2 , x3 ] f [ x1 , x2 , x3 ] f [ x2 , x3 ] f [ x1 , x2 , x3 , x4 ] f [ x2 , x3 , x4 ] f [ x3 , x4 ] f [ x0 , x1 ,?, x4 ]

x2 f ( x2 ) x3 f ( x3 ) x4 f ( x 4 )

规定函数值为零阶均差
252

5.3.1 均差概念及计算
均差表的数据构成一个矩阵F: F00=f(x0) F10=f(x1),F11=f [x0,x1 ] F20=f(x2),F21=f [x1,x2],F22=f [x0,x1,x2] F30=f(x3),F31=f[x2,x3],F32=f[x1,x2,x3], F =f [ x ,x ,x ,x ] 33 0 1 2 3 一般有 Fi, j=f [xi-j , xi-j+1 ,…,xi-1 ,xi ] Fi, j-1=f [xi-j+1 ,…,xi ]

Fij

Fi-1, j-1=f [xi-j ,,…,xi-1 ]

253

5.3.1 均差概念及计算
计算机上计算均差表的公式
i ? 0,1,...,n ? Fi 0 ? f ( x i ) , ? ? Fi j ? ( Fi j ?1 ? Fi ?1 j ?1 ) /( x i ? x i ? j ),

j ? 1,2,...,n.

i ? j , j ? 1,...,n

254

5.3.1 均差概念及计算
例5-4 已知函数y=f(x)的观测数据如表,试构造 均差表,并求f [2,4,5]及f [2,4,5,6]的值。 解 n=4, 构造均差表
xi 0 2 4 f(xi ) 1阶 1 5 9 2 2阶 3阶 4阶

x

0

2

4

5

6

f(x) 1

5

9

- 4 13

5
6

-4
13

5?1 ?2 2?0 2 9 ? 50? 2 2 ? 2 ? 0 4 0? 2 ? 5 ? 0 ?? 42 ?9 4 ?? 13 ? ?1 ? ?13 ?5 ?5 -13 5 ?-5 -1 ?0 4 5?2 ? ( ?4 ) 17 ? ( ?13)15 ? ( ?5)5 ? ( ?1) 17 1315 1 ? 15 ? 5 ? 1 ?5 17 6?5 6?4 6?2 6?0

f[2,4,5]= -5 f[2,4,5,6]=5
255

5.3.1 均差概念及计算(回顾)
目的:构造具有如下形式的插值多项式:

P n ( x) = a0 + a1 ( x - x0 ) + a2 ( x - x0 )( x - x1 ) + L + an ( x - x0 )( x - x1 ) L ( x - xn- 1 )
它满足递推性:
P n ( x) = P n- 1 ( x) + an ( x - x0 )( x - x1 ) L ( x - xn- 1 )

256

5.3.1 均差概念及计算(回顾)
零阶均差定义为函数值本身,即 f [ xi ] ? f ( xi ) 一阶均差

x j ? xi (一阶差商 ) f [ x j , xk ] ? f [ xi , x j ] 二阶均差 f [ xi , x j , xk ] ? xk ? xi k阶均差
f [ xi ?k ?1 , ?, xi ?1 , xi ] ? f [ xi ?k , xi ?k ?1 , ?, xi ?1 ] f [ x i ? k , x i ? k ?1 , ? , x i ] ? xi ? xi ?k
257

f [ xi , x j ] ?

f ( x j ) ? f ( xi )

5.3.1 均差概念及计算 粗线框出的部分在计算机
均差表1
xi f(xk) 1阶 2阶

上可存入二维数组
3阶
f [ x 2 , x 3 , x4 ] ?

4阶
f [ x 3 , x4 ] ? f [ x 2 , x 3 ] x4 ? x 2

x0 f(x0)
x1 f(x1) f [x0,x1 ] x2 f(x2) f [x1,x2] f [x0,x1,x2]

f [ x2 , x3 , x4 ] ? f [ x1 , x2 , x3 ] f [ x1 , x2 , x3 , x4 ] ? x4 ? x1

x3 f(x3) f [x2,x3 ] f [x1,x2,x3 ] f [x0,x1,x2,x3 ]
x4 f(x4) f [x3,x4] f [x2,x3,x4 ] f [x1,x2,x3,x4 ] f [x0,x1,x2,x3,x4 ]

┊ ┊





┊ (回顾) ┊
258

5.3.1 均差概念及计算(回顾)
计算机上计算均差表的公式
i ? 0,1,...,n ? Fi 0 ? f ( x i ) , ? ? Fi j ? ( Fi j ?1 ? Fi ?1 j ?1 ) /( x i ? x i ? j ),

Fij
j ? 1,2,...,n. i ? j , j ? 1,...,n

259

5.3.2牛顿均差型插值多项式
f ( x1 ) - f ( x0 ) ( x - x0 ) 根据线性插值的点斜式 P1 ( x) = f ( x0 ) + x1 - x0

可得牛顿均差型线性插值多项式:

P1 ( x) ? f ( x0 ) ? f [ x0 , x1 ]( x ? x0 )

P2 ( x) = P 1 ( x) + a2 ( x - x0 )( x - x1 )
= f ( x0 ) + f [ x0 , x1 ]( x2 - x0 ) + a2 ( x2 - x0 )( x2 - x1 )

由 f ( x2 ) = P2 ( x2 ) 可得 a2 = ( f [ x0 , x2 ] - f [ x0 , x1 ]) /( x2 - x1 ) = f [ x0 , x1, x2 ] 牛顿均差型二次插值多项式:

P2 ( x) ? f ( x0 ) ? f [ x0 , x1 ]( x ? x0 ) ? f [ x0 , x1 , x2 ]( x ? x0 )( x ? x1 ) 260

5.3.2牛顿均差型插值多项式
牛顿均差型n次插值多项式(n次牛顿均差插值多项 式) Pn ( x) ? f ( x0 ) ? f [ x0 , x1 ]( x ? x0 ) ? f [ x0 , x1 , x2 ]( x ? x0 )( x ? x1 )
? ?? f [ x0 , x1 ,?, xn ]( x ? x0 )( x ? x1 )?( x ? xn?1 )

余项公式的牛顿形式:
Rn ( x) ? f [ x, x0 , x1 ,?, xn ]( x ? x0 )( x ? x1 )?( x ? xn )

计算牛顿均差插值多项式的步骤: (1)作均差表 (2)根据公式计算牛顿型插值多项式(表中 对角线上各均差值就是Pn(x)的各项系数)。

261

5.3.2牛顿均差型插值多项式
例5-5 已知函数y=f(x)的观测数据如例5-4,试用 全部基点构造牛顿均差插值多项式,并用二次插 值求f(3)的近似值。 解 用全部基点时,n=4,先作均差表,见例5-4。
x 0 2 4 5 6
f(x) 1 5 9 - 4 13
xi 0 2 4 f(xi ) 1阶 1 5 9 2 2 -13 17 0 -5 15 2阶 3阶 4阶

5
6

-4
13

-1
5 1
262

5.3.2牛顿均差型插值多项式
P4(x)=f(0)+f[0,2](x- 0)+f[0,2,4](x- 0)(x- 2)+ f[0,2,4,5](x- 0)(x- 2)(x- 4)+f[0,2,4,5,6](x- 0)(x- 2)(x- 4)(x- 5) =1+2x- x(x- 2)(x- 4)+x (x- 2)(x- 4)(x- 5)

用二次插值求 ff(x (3) 时 n =2 , x=3 ,作内插取 xi ) 1 阶 2 阶 3 阶 4阶 i
0 1 =4, x =5 x0=2, x 1 2

f(3)≈P2(3)=f(2)+ 2 0 – 2)+ f[2,4,5](3- 2)(3- 4) 4 9 f[2,4](3
=7- 5(3- 2)(3- 4) = 12 17 15 5 1 6 13 ?若插值基点等距分布,牛顿型插值多项式还可以 利用差分得到简化。 263
5 -4 -13 -5

2

5

2

-1

5.3.2牛顿均差型插值多项式
?差分 设基点等距, xi=x0+ ih, i=0,1,…,n. h称为步长。 一阶差分: xi+1 = xi + h 一阶向前差分、一阶前差 D f ( xi ) = f ( xi+ 1 ) - f ( xi ) 一阶向后差分、一阶后差 ? f ( xi ) f ( xi ) - f ( xi- 1 ) 二阶差分: 二阶前差 D 2 f ( xi ) = D (D f ( xi )) = D f ( xi+ 1 ) - D f ( xi ) 二阶后差 ? 2 f ( xi )
蜒 ( f ( xi )) = ? f ( xi ) f ( xi- 1 )
264

5.3.2牛顿均差型插值多项式
k阶前差和k阶后差
D k f ( xi ) = D (D k- 1 f ( xi )) = D k- 1 f ( xi+ 1 ) - D k- 1 f ( xi )
? k f ( xi ) 蜒 (
k- 1

f ( xi )) = ?

k- 1

f ( xi )

k- 1

f ( xi- 1 )

零阶差分规定为 D 0 f ( xi ) = ? 0 f ( xi )
D k f ( xi ) ,
k

fi = f ( xi )
k

f ( xi )

简记为

D fi ,

k

fi

265

5.3.2牛顿均差型插值多项式
差分表
xk f k 一n 差分 x0 f 0 x1 f1 二n 差分 三n 差分 四n 差分

?f 0
?f 1

?f 1 ?f 2

?2 f0

? 2 f2
?2 f 3
?2 f 4

?3 f0
3

x2 f 2

?f 2
x3 f3 x4 f 4

?2 f1 ? f2
2

?3 f 3

? f1
?3 f 4

?4 f0
?4 f4

?f 3
?f 4

?f 3

266

5.3.2牛顿均差型插值多项式
例5-6 设f(x)=x2+ x ,取步长 h=0.5,计算D 2 f (0) 解 令 x0=0, 则 x1=0.5, x2=1 △f(0)=f(0.5)- f(0) =0.75 , △f(0.5)=f(1)- f(0.5) =1.25 △2f(0)=△(△f(0))=△f(0.5)- △f(0) =1.25 – 0.75 = 0.5
? f (1) f (1) - f (0.5) = 1.25
,
2

f (1)

? f (0.5)

f (0.5) - f (0) = 0.75
1.25 - 0.75 = 0.5
267

? 2 f (1)

蜒 ( f (1)) = ? f (1) ? f (0.5)

5.3.2牛顿均差型插值多项式
用数学归纳法可以证明前差、后差、均差之间的 下述关系:
(1) ? f i ? ? f i ? k
k k

1 1 k k ( 2) f [ x0 , x1 ,? , xk ] ? ? f ? ? fk 0 k k k! h k! h

利用差分可以简化等距基点的牛顿型插值公式, 称为前差(或后差)插值公式 。

268

D 3 yn = D (D 2 yn ) = D 2 yn+ 1 - D 2 yn = D (D yn+ 1 ) - D (D yn ) = D yn+ 2 - D yn+ 1 - (D yn+ 1 - D yn ) = yn+ 3 - yn+ 2 - 2( yn+ 2 - yn+ 1 ) + yn+ 1 - yn = yn+ 3 - 3( yn+ 2 - yn+ 1 ) - yn =2 =2
n+ 3

- 3(2

n+ 2

- 2 )- 2
n

n+ 1

n

= (8 - 12 + 6 - 1) 2
n

269

5.4 三次Hermite插值
?5.4.1
?5.4.2

两点三次Hermite插值
两点三次Hermite插值的余项

?5.4.3

分段两点三次Hermite插值

270

5.4 Hermite插值
? Newton插值和Lagrange插值要求插值多项式满 足条件 P( xi ) = f ( xi ) = yi i = 0,1,L , n --------(1) ?Newton插值和Lagrange插值虽然构造比较简单, 但不能很好的逼近函数f(x)。

?Newton插值和Lagrange插值的高次插值稳定性差, 对于计算过程的舍入误差十分敏感。

271

5.4 Hermite插值
?f(x)在区间[a,b]上具有一阶导数,设P(x)为f(x) 在区间[a,b]上的插值函数。若要求P(x)在[a,b] 上具有一阶导数(一阶光滑度)即P(x)在插值基点 x0,x1,……xn处必须满足:
P( xi ) = f ( xi ) = yi
Pⅱ ( xi ) = f ( xi ) = yi

i = 0,1,L , n

--------(1) --------(2)

i = 0,1,L , n

满足(1)和(2)式的插值问题为Hermite插值。 称满足(1)和(2)式的插值多项式P(x)为Hermite插 值多项式,记为Hk(x) , k为多项式次数。

272

5.4.1 两点三次Hermite插值
考虑只有两个插值基点的插值问题: 设f(x)在插值基点x0,x1处的函数值为y0,y1 在x0,x1处的一阶导数值为y0',y1' 两个插值基点x0,x1最高可以用3次Hermite多项式 H3(x)作为插值函数。 H3(x)=a0+a1x+a2x2+a3x3

H3(x)应满足插值条件: H3 ( x0 ) = y0
ⅱ H3 ( x0 ) = y0

H3 ( x1 ) = y1
ⅱ H3 ( x1 ) = y1
273

5.4.1 两点三次Hermite插值
H3(x)=a0+a1x+a2x2+a3x3,满足插值条件:
H3 ( x0 ) = y0
ⅱ H3 ( x0 ) = y0
ì H 3 ( x0 ) = y0 ? ? ? ? H 3 ( x1 ) = y1 ? ? í ? ? ? H 3 ( x0 ) = y0 ? ? ? ? ?( x1 ) = y1? H ? 3 ? ?

H3 ( x1 ) = y1
ⅱ H3 ( x1 ) = y1
y0 y1 ? y0 y1?
274

证明满足以上插值条件的H3(x)存在且唯一。
2 3 ì ? a + a x + a x + a x = 0 1 0 2 0 3 0 ? ? 2 3 ? a + a x + a x + a x ? 0 1 1 2 1 3 1 = ? í 2 ? a + 2 a x + 3 a x ? 1 2 0 3 0 = ? ? 2 ? a + 2 a x + 3 a x 1 2 1 3 1 = ? ?

5.4.1 两点三次Hermite插值
系数矩阵:

275

5.2.3 n次拉格朗日插值
作基点xi的n次插值基函数 ( i=0,1,…,n):
( x ? x0 )( x ? x1 )?( x ? xi ?1 )( x ? xi ?1 )?( x ? xn ) li ( x ) ? ( xi ? x0 )( xi ? x1 )?( xi ? xi ?1 )( xi ? xi ?1 )?( xi ? xn ) ??
n

x ? xj xi ? x j

j ?0 j?i

i ? 0,1,? , n

基函数具有性质: li (xi ) = 1 , li (x j ) = 0 n次拉格朗日插值多项式

(j ? i),
n

i、j
n

0,1,L , n
x ? xj xi ? x j
278

Pn ( x) ? y0l0 ( x) ? y1l1 ( x) ? ? ? ynln ( x) ? ? yi ?
i ?0 j ?0 j?i

5.4.1 两点三次Hermite插值
希望插值系数与Lagrange插值一样简单!
? H3(x)应使用四个插值基函数表示,设H3(x)的插 值基函数为 hi ( x), i = 0,1, 2,3

H3 ( x) = a0h0 ( x) + a1h1 ( x) + a2h2 ( x) + a3h3 ( x)
重新假设

ⅱ H3 ( x) = y0a 0 ( x) + y1a 1 ( x) + y0 b 0 ( x) + y1b1 ( x)

ⅱ ⅱ H3 ( x) = y0a 0 ( x) + y1a 1ⅱ ( x) + y0b0 ( x) + y1b1 ( x) 279

5.4.1 两点三次Hermite插值
其中
?0 ( x0 ) ? 1 ?1 ( x0 ) ? 0
? 0 ( x0 ) ? 0 ? 1 ( x0 ) ? 0

?0 ( x1 ) ? 0 ?1 ( x1 ) ? 1
? 0 ( x1 ) ? 0 ? 1 ( x1 ) ? 0

? ( x0 ) ? 0 ?0 ? ( x0 ) ? 0 ?1 ? ( x0 ) ? 1 ?0 ? ( x0 ) ? 0 ?1

? ( x1 ) ? 0 ?0 ? ( x1 ) ? 0 ?1 ? ( x1 ) ? 0 ?0 ? ( x1 ) ? 1 ?1

可知x1是α0(x)的二重零点,即可假设

a 0 ( x) = ( x - x1 ) (ax + b)

a 0 ( x0 ) = 1

2

?( x0 ) = 0 a0

280

5.4.1 两点三次Hermite插值
可得
2 a= ( x0 - x1 )3 2 x0 1 b= + 2 ( x0 - x1 ) ( x0 - x1 )3

a 0 ( x) = ( x - x1 )2 (ax + b)

? 2x 2 ?= ( x - x1 ) ? 3 ? è ( x0 - x1 )

2 x0 ? 1 ÷ ÷ + + 2 3÷ ( x0 - x1 ) ( x0 - x1 ) ÷ ?

( x - x1 )2 = ( x0 - x1 )2

? ? 2 x0 2 x ? ÷ 1 + ? ÷ ? ÷ è x0 - x1 x0 - x1 ÷ ?

2 骣 骣 x - x0 ÷ ? x - x1 ÷ ? ÷ ÷ ? =? 1 + 2 ÷ ÷ ? ? ÷ x x x1 - x0 桫0 1 ÷ 桫

Lagrange 插值基函数
281

2 = (1+ 2l1 ( x)) l0 ( x)

5.4.1 两点三次Hermite插值


类似可得:

2 骣 骣 x - x1 ÷ x x 2 ? ÷ 0 ? a 0 ( x) = (1+ 2l1 ( x)) l0 ( x) = ?1 + 2 ÷ ÷ ? ÷ ÷ ? ? ÷ x x x1 - x0 桫0 1 ÷ 桫 2

骣 x - x0 ÷ x - x1 ÷ 骣 ? ? ÷ ÷ 1 + 2 ? a 1 ( x) = (1+ 2l0 ( x)) l ( x) = ? ÷ ÷ ? ? ÷ x x x0 - x1 桫1 0 ÷ 桫
2 1 2 ( x) = (x - x0 ) b 0 ( x) = ( x - x0 ) l0

骣 x - x1 ÷ ? ÷ ? ÷ ? ÷ x x 桫0 1 骣 x - x0 ÷ ? ÷ ? ÷ ? ÷ x x 桫1 0

2

2

b1 ( x) = ( x - x1 ) l ( x) = (x - x1 )
2 1

将以上结果代入
ⅱ H3 ( x) = y0a 0 ( x) + y1a 1 ( x) + y0 b0 ( x) + y1b1 ( x)
282

5.4.1 两点三次Hermite插值
得到两点的三次Hermite插值公式

ⅱ H3 ( x) = y0a 0 ( x) + y1a 1 ( x) + y0 b 0 ( x) + y1b1 ( x)
2 = y0 (1+ 2l1 ( x)) l0 ( x) + y1 (1+ 2l0 ( x)) l12 ( x)

2 2 ? ? + y0 ( x - x0 ) l0 ( x) + y1 ( x - x1 ) l1 ( x)

? ? x ? x0 ? ? x ? x1 ? x ? x1 ? ? x ? x0 ? ? y0 ? ?x ?x ? ? ? y1 ? ?1 ? 2 x ? x ? ?? ?1 ? 2 x ? x ? ?? ?x ?x ? ? 1 ? 0 ? 1 0 ?? 0 0 1 ?? 1 ? ?
? x ? x0 ? ? x ? x1 ? ? ? ?x ? x1 ? ? ? ?x ? x0 ? ? ? y1 ? y0 ? ? ? ?x ?x ? x ? x 0 ? ? 1 1 ? ? 0
2 2

2

2

283

5.4.2 两点三次Hermite插值的余项
两点三次Hermite插值的误差为
R3 ( x) = f ( x) - H3 ( x)

R3 ( xi ) = f ( xi ) - H3 ( xi ) = 0

i = 0,1

ⅱ R3 ( xi ) = f ( xi ) - H3 ( xi ) = 0

x0,x1均为R3(x)的二重零点,因此可设
R3 ( x) = K ( x)( x - x0 )2 ( x - x1 )2

其中K(x)待定

284

5.4.2 两点三次Hermite插值的余项
构造辅助函数
j (t ) = f (t ) - H3 (t ) - K ( x)(t - x0 ) (t - x1 )
2 2

均是 二重根
i = 0,1

j ( xi ) = f ( xi ) - H3 ( xi ) - K ( x)( xi - x0 )2 ( xi - x1 )2 = 0

j ( x) = f ( x) - H3 ( x) - K ( x)( x - x0 )2 ( x - x1 )2 = 0

因此j ( x) 至少有3个零点:x、x0 (二重) 、x1(二重) 连续使用4次Rolle定理,可得至少存在一点

x ? [ x0 , x1 ] 使得

j

(4)

(x) = 0

285

5.4.2 两点三次Hermite插值的余项


j (t )(4) = f (t ) - (4) H3 (t ) - K ( x)(t - x0 ) (t - x1 )
j (x) = f (x) - 4! K ( x) = 0
f (4) (x ) K ( x) = 4!

2

2

所以,两点三次Hermite插值的余项为
f (4) (x ) R3 ( x) = ( x - x0 ) 2 ( x - x1 ) 2 4!
a# x b

x0 #x

x1

(4) (4) max f ( x ) = M 代替 f (x ) 以上分析都能成立吗?

当f (4)(x)在[x0,x1]上存在时,上述余项公式成立!
286

5.4.2 两点三次Hermite插值的余项
例5-7 f(x)在插值基点1,2处的函数值为f(1)=2, f(2)=3, f(x)在插值基点1,2处的导数值为 f '(1)=0,f '(2)= -1,求f(x)的两点三次插值多项式 及f(x)在x=1.5,1.7处的函数值。
ⅱ = 0, y1 = - 1 解: x0 = 1, x1 = 2 y0 = 2, y1 = 3 y0

ⅱ H3 ( x) = y0a 0 ( x) + y1a 1 ( x) + y0 b 0 ( x) + y1b1 ( x)
2 骣 骣 骣 骣 x - x0 ÷ ? x - x1 ÷ x - x1 ÷ ? x - x0 ÷ ? ? ÷ ÷ ÷ ÷ ? ? = y0 ? 1 + 2 + y 1 + 2 ? ÷ ÷ 1 ÷ ÷ ? ? ? ? ÷ ÷ ÷ ÷ x x x x 桫 桫 x x x x 桫 桫 0 1 1 0 1 0 0 1 2

骣 x - x1 ÷ ? ? ÷ + y0 (x - x0 ) ? ÷ ? ÷ x x 桫0 1

2

骣 x - x0 ÷ ? ÷ + y1? (x - x1 ) ? ? ÷ x1 - x0 ÷ 桫

2

287

5.4.2 两点三次Hermite插值的余项
H3 ( x) = 2(1+ 2( x - 1))( x - 2) + 3(1- 2( x - 2)) ( x - 1)
2 2

- ( x - 2) ( x - 1)

2

= - 3x3 + 13x 2 - 17 x + 9 f (1.5) ? H3 (1.5) = 2.625 f (1.7) ? H3 (1.7) = 2.931

作为多项式插值,三次已是较高的次数,次数再 高就有可能发生Runge现象,因此,对有n+1节点 的插值问题,可以使用分段两点三次Hermite插值 288

5.4.3 分段两点三次Hermite插值

定 义 : 设 函 数 f(x) 在 n+1 个 插 值 基 点 a≤x0<x1<……<xn≤b处的函数值和导数值为f(xi )=yi , f '(xi )=mi ,i=0,1,2,…,n,满足如下条件的函数称之 为分段三次Hermite插值函数:

(1)H(x)在每个子区间[xk,xk+1]上是不超过三次
的多项式。

(2) H(xi )=yi , H '(xi )= mi i=0,1,2,…,n
(3) H '(x)在区间[a , b]上连续.
289

5.5 三次样条插值

?5.5.1

三次样条插值函数的概念

?5.5.2

三次样条插值函数的求法

293

5.5.1 三次样条插值函数的概念
样条插值的思想:逐段选取适当的低次多项 式,按一定的光滑性要求连接起来构成插值 函数。

294

5.5.1 三次样条插值函数的概念

定义:给定区间[a,b]上n+1个点a=x0<x1<……<xn=b, 以及相应的函数值 yi=f(xi), i=0,1,…,n 。如果函数 S(x)满足: (1)在每个子区间 [xk , xk+1](k=0,1,…,n-1)上, S(x) 是不超过三次的多项式,且S(xi )=yi ,i=0,1,…,n;

(2) S ( x), S ⅱ ( x), S ( x) 在[a, b]上连续;
则称S(x)是f(x)在基点x0 , x1 , x2 ,…, xn上的三次样条 插值函数. 称xoy平面上的点(xi , yi )(i=0,1,…,n)为样点。
295

5.5.1 三次样条插值函数的概念
例 5-8 给 定 区 间 [0, 3] 上 3 个 点 的 函 数 值 f(0)=0, f(1)=2, f(3)=4,试求数 a,b,c,d,使函数 S(x)为给定点 上的三次样条插值函数。 解:
2 ì ? x + x + d, ? S ( x) = í 3 2 ? ax + bx + cx + 1, ? ?

0 #x 1 #x

1 3

S1 ( x) = x + x + d

2

S2 ( x) = ax3 + bx2 + cx + 1
S1ⅱ ( x) = 2,
ⅱ S2 ( x) = 6ax + 2b,
296

则 S ? ( x) = 2x + 1, 1
? ( x) = 3ax2 + 2bx + c, S2

5.5.1 三次样条插值函数的概念
根据定义,由
由 由 由 由
S (0) = S1 (0) = f (0)

得 d=0,故

S1 ( x) = x2 + x

S1 (1) = S2 (1) = f (1)

S2 (3) = f (3)

得 a+ b+ c = 1 得 27a + 9b + 3c = 3 得 得
3a + 2b + c = 3

S1ⅱ (1) = S2 (1)

ⅱ S1ⅱ (1) = S2 (1)

6a + 2b = 2

求得

a= - 1 , b= 4 , c= - 2 , d = 0 2 ì ? x + x, 0 #x 1 ? S ( x) = í 3 2 ? x + 4 x - 2 x + 1, 1 #x 3 ? ?

297

5.5.1 三次样条插值函数的概念
给定n+1个样点(xi , yi )(i=0,1,…,n) ,确定一个 三次样条插值函数需要4n个独立条件。在定义中, 已指定了4n-2个条件,即
ì S ( xi ) = yi i = 0, 1, 2,...n ? ? ? ? S- ( xi ) = S+ ( xi ) = yi , i = 1, 2,...n - 1 ? ? í ? S-ⅱ ( xi ) = S+ ( xi ), i = 1, 2,...n - 1 ? ? ? ⅱ( xi ), ? i = 1, 2,...n - 1 ? ? S-ⅱ( xi ) = S+

所以,一般需补充指定2个边界条件。
298

5.5.2 三次样条插值函数的求法
常用的边界条件有三种: (1)第一种边界条件:给定端点处的一阶导数,即
Sⅱ ( x0 ) = y0 , S ⅱ ( xn ) = yn;

(2)第二种边界条件:给定端点处的二阶导数m0 ,mn。 特别取 m0 = mn=0 时称为自然边界条件,求得的 S(x)称为自然样条函数。 (3)第三种边界条件:当 y=f(x)是以 b- a为周期的周 期函数时,端点要满足 S ? ( x ) ? S ? ( x ), m ? m 。 304
0 n 0 n

5.6 B样条插值
?5.6.1
?5.6.2

B样条函数
B样条曲线

?5.6.3

自由曲线设计

311

5.6.1 B样条函数
在我们工程中应用的拟合曲线,一般地说可 以分为两种类型:一种是最终生成的曲线通过所 有的给定型值点,比如抛物样条曲线和三次参数 样条曲线等,这样的曲线适用于插值放样;另一 种曲线是,它的最终结果并不一定通过给定的型 值点,而只是比较好地接近这些点,这类曲线 (或曲面)比较适合于外形设计。

312

5.6.1 B样条函数
因为在外形设计中 ( 比如汽车、船舶 ) ,初始 给出的数据点往往并不精确;并且有的地方在外 观上考虑是主要的 ,因为不是功能的要求,所以 为了美观而宁可放弃个别数据点。因此不须最终 生成的曲线都通过这些数据点。另一方面,考虑 到在进行外形设计时应易于实时局部修改,反映 直观,以便于设计者交互操作。第一类曲线在这 方面就不能适应。
313

5.6.1 B样条函数
到 了 70 年 代 , 法 国 雷 诺 汽 车 公 司 的 工 程 师 贝 齐 尔 ( Bezier )创造出一种适用于几何体外形设计的新的曲 线表示法。这种方法的优越性在于:对于在平面上随手 勾画出的一个多边形(称为特征多边形),只要把其顶 点坐标输入计算机,经过不到一秒钟的计算,绘图机就 会自动画出同这个多边形很相像、又十分光滑的一条曲 线。这种方法被人们称为贝齐尔 (Bezier) 方法 ( 以下统称 为Bezier方法)。

314

5.6.1 B样条函数
贝齐尔曲线的形状是通过一组多边折线(也称 为贝齐尔控制多边形)的各顶点惟一地定义出来 的。在该多边折线的各顶点中,只有第一点和最 后一点在曲线上,其余的顶点则用来定义曲线的 形状。列举一些 Bezier 多边折线和相应的 Bezier 曲线的形状关系如图5.6-1所示。

315

5.6.1 B样条函数

图5.6-1 Bezier

曲线
316

5.6.1 B样条函数

317

5.6.1 B样条函数

称为n次Bernstein多项式的基函数,Bezier曲线就 是以此为基础构造出来的。 318

5.6.1 B样条函数

图5.6-2 2到8次Bezier曲线图例

319

5.6.1 B样条函数

为了定义 B样条曲线,首先给出 n次截幂函数和 n 阶B样条函数的定义。我们称 为n次截幂函数, 即

称Mn(x)为n阶B样条函数,即

320

5.6.1 B样条函数

图5.6-3

B样条函数

321

5.6.1 B样条函数
?B样条函数具有下列的重要性质:
?(1)Mn(x)是分段n-1次多项式。当n为偶数时

具有整数节点 xk=-n/2+k ,当 n 为奇数时具有半 整数节点:xk=-n/2+k。比如:
?n=2, ?n=3,

k=0,1,2, x0=-1, x1=0, x2=1 k=0,1,2,3,x0=-1.5, x1=-0.5, x2=0.5,

x3=1.5
322

5.6.1 B样条函数
?B样条函数具有下列的重要性质:
?( 2 )在整个实数轴上, Mn(x) 具有 n-2 次连续

导数。
?(3)Mn(x)是偶函数。 ?( 4 ) Mn(x) 具有非负性质,即当 |x| < n/2 时,

Mn(x)>0;当|x|>n/2时,Mn(x)= 0;且0≤Mn(x)≤1。

323

5.6.2 B样条曲线
给定一组初始型值点Pi(i=0,1,…,n),将它们按次 序连接为折线段,称为控制多边形。我们称

为m次B样条曲线。它是对参数t具有m-1阶连续导 数的分段m次多项式。
324

5.6.2 B样条曲线

图5.6-4 一次B样条曲线S1(t)
325

5.6.2 B样条曲线

图5.6-5 二次B样条曲线S2
326

5.6.3 自由曲线设计
(1)设计的曲线要求出现直线段。

(2)设计的曲线要求出现尖点。只需在出现尖 点处将型值点的序号重复编号。
下面给出在计算机上绘制二次B样条曲线的步骤: (1) 首先要确定全部的型值点Pi (i=0,1,…,n)。 由于Pi (i=0,1,…,n)是平面上的点,可通过定 义二维数组向Pi (i=0,1,…,n)赋坐标值。
327

5.6.3 自由曲线设计

(2) 用直线连接型值点Pi (i=0,1,…,n),画出控制多 边形。 (3) 整条曲线需分n-1段绘制。对于第i (i=0,1,…, n-2)段,其型值点为Pi,Pi+1,Pi+2。给出区间[0, 1]上的加密点数m,则u=j/m,计算出各加密点的 坐标B2(Pi,Pi+1,Pi+2,j/m)(j=0,1,…,m),将 这m+1个加密点用直线连接,绘制出本段B样条曲 (4) 所有的分段曲线全部绘出,即完成了整条曲线 328 的绘制。

5.7 曲线拟合的最小二乘法

?5.7.1

曲线拟合的最小二乘法

?5.7.2

代数多项式拟合

329

5.7.1 曲线拟合的最小二乘法
曲线拟合的问题: 设函数y=f(x)在n个互异点的观测数据为 xi yi x1 y1 x2 ….. xn y2 ….. yn

求一个简单的近似函数φ(x),使之 “最好” 地逼近f(x),而不必满足插值原则。称函数 y=φ(x)为经验公式或拟合曲线。
330

5.7.1 曲线拟合的最小二乘法
例子: (注意它与插值法的不同)

331

5.7.1 曲线拟合的最小二乘法
通常选择函数类型的做法:描出散点图,再根 据专业知识和经验来选择φ(x)的类型。

? ( x) ? c1?1 ( x) ? c2? 2 ( x) ? ?? cm? m ( x),(m ? n)
c1 , c2 ,?, cm 是待定常数,

j i ( x)(i = 1, 2,L , m)
为某类函数中线性无关的简单函数。
332

5.7.1 曲线拟合的最小二乘法
如取 ? j ( x) ? x
j ?1

就得到代数多项式
2 m ?1

? ( x) ? c1 ? c2 x ? c3 x ? ?? cm x
近似曲线在各点的偏差的平方和是
n

2 2 ( ? ( x ) ? y ) ? ( c ? ( x ) ? c ? ( x ) ? ? ? c ? ( x ) ? y ) ? i i ? 1 1 i 2 2 i m m i i i ?1

n

i ?1
n

? ? (ai 1c1 ? ai 2c2 ? ? ? aim cm ? yi )2
i ?1

ai j ? ? j ( x i )

333

5.7.1 曲线拟合的最小二乘法
线性最小二乘法: 求近似函数 ? ( x) ? c1?1 ( x) ? c2? 2 ( x) ? ?? cm? m ( x),(m ? n) 使其在指定n个点的偏差平方和达最小。 问题转化为求待定系数
n i ?1

c1 , c2 ,L , cm

使 ? (ai1c1 ? ai 2c2 ? ? ? aimcm ? yi )2 达到极小。 其中 ai j ? ? j ( xi )
334

5.7.1 曲线拟合的最小二乘法
考虑超定方程组

ì a11c1 + a12 c2 + L + a1m cm = y1 ? ? ? ? LLLLLL í ? ? ? ? ? an1c1 + an 2 c2 + L + anm cm = yn


(n>m)

Ac = y
T

A = [ai j ]n? m

c = [c1, c2 ,L , cm ]

y = [ y1 , y2 ,L , yn ]

T
335

5.7.1 曲线拟合的最小二乘法
如果有向量c使得
Ac ? y 2 ? min
2
2

? (a
i ?1

n

i1 1

c ? ai 2c2 ? ? ? aimcm ? yi )

达到最小,则称c为超定方程组的最小二乘解。 定理 5-3 超定方程组Ac = y 存在最小二乘解, 且即为方程组ATAc = ATy 的解。当A的列向量线 性无关时ATA非奇异,这时有唯一的解。
336

5.7.1 曲线拟合的最小二乘法
称方程组ATAc = ATy 为方程组Ac = y 的正则方程 组、正规方程组、法方程组. 曲线拟合的最小二乘法可以看成求下述超定方程 组的最小二乘解的问题:
ì c1j 1 ( x1 ) + c2j 2 ( x1 ) + L + cmj m ( x1 ) = y1 ? ? ? ? ? c1j 1 ( x2 ) + c2j 2 ( x2 ) + L + cmj m ( x2 ) = y2 í ? LLLLLL ? ? ? ? ? ? c1j 1 ( xn ) + c2j 2 ( xn ) + L + cmj m ( xn ) = yn

简写为

? ( xi ) ? yi ,

i ? 1,2,...,n

337

5.7.1 曲线拟合的最小二乘法
一般计算步骤为: (1)假设经验公式
? ( x) ? c1?1 ( x) ? c2?2 ( x) ?
? cm?m ( x)

计算 A = [aij ]n? m ,其中

aij = j j ( xi ) , i = 1,2,L , n, j = 1,2,L , m
(2)计算ATA, ATy ,形成法方程组ATAc = ATy

(3)求解法方程组,输出 c1 ,c2 ,…,cm ,构成

? ( x) ? c1?1 ( x) ? c2?2 ( x) ?

? cm?m ( x)
338

5.7.1 曲线拟合的最小二乘法
例 已知观测数据(1,–5),(2,0),(4,5),(5,6) ,试用最
b 小二乘法求形如 j ( x) = ax + x 的经验公式。 解 j ( x) = aj 1 ( x) + bj 2 ( x) j 1 ( x) = x , j 2 ( x) = 1/ x
ì j ? ? 作超定方程组 ? ? ?j í ? j ? ? ? ? ?j
轾 1 1 犏 犏 2 0.5 犏 A= 犏 4 0.25 犏 犏 5 0.2 臌

(1) = - 5 (2) = 0 (4) = 5 (5) = 6



ì a + b= - 5 ? ? ? ? ? 2a + 0.5 b = 0 í ? 4a + 0.25b = 5 ? ? ? ? ? 5a + 0.2b = 6

轾 - 5 犏 犏 0 y= 犏 犏 5 犏 犏 6 臌

339

5.7.1 曲线拟合的最小二乘法
轾 46 4 A A= 犏 犏 4 1.3525 臌
T

轾45 A y= 犏 犏 - 2.55 臌
T

ì 46a + 4b = 45 ? ? 法方程组为í ? ? ? 4a + 1.3525b = - 2.55

求得 a=1.537650114 b= -6.432976311
6.432976 所求经验公式为 j ( x) = 1.537650 x x
340

5.7.2 代数多项式拟合
1、直线拟合.? ( x ) ? c1 ? c2 x m = 2, j ( x) = 1, j ( x) = x 1 2 作超定方程组
ì c1 + c2 x1 = y1 ? ? ? ? c1 + c2 x2 = y2 ? í ? LLLLLL ? ? ? ? ? ? c1 + c2 xn = yn

轾 1 犏 犏 1 犏 A= 犏 M 犏 犏 1 臌

x1 x2 M xn

记号?指 对i从1到 n取和

轾 1 1 A A= 犏 犏 x1 x2 臌
T

轾 1 犏 ... 1 犏 1 犏 ... xn 犏 M 犏 犏 1 臌

x1 x2 M xn =

n

?

xi

? ?

xi xi2
341

5.7.2 代数多项式拟合
轾 1 1 A y= 犏 犏 x1 x2 臌
T

轾 y1 犏 ... 1 犏 y2 犏 = ... xn 犏 M 犏 犏 yn 臌

?

?

yi

xi yi

法方程组为
? ? c1n ? c2 ? xi ? ? yi ? 2 c x ? c x ? 2 ? i ? ? x i yi ? 1? i
342

5.7.2 代数多项式拟合
2、二次拟合、抛物拟合 ? ( x) ? c1 ? c2 x ? c3 x 作超定方程组 轾 ì ?c + c x + c x = y 1 x x2
1 ? 1 ? ? 2 ? ? c1 + c2 x2 + c3 x 2 = y2 í ? LLLLLL ? ? ? 2 ? c + c x + c x 2 n 3 n = yn ? ? 1 2 1 3 2

2

轾n 犏 T A A= 犏 邋xi 犏 犏 x2 邋i 犏 臌

邋x
x x

i 2 i 3 i

xi2 xi3 xi4

1 犏 1 2 犏 1 x x 有 A= 犏 2 2 犏 M M M 犏 犏 2 1 x x 犏 n 臌 n 轾 ? yi 犏 AT y = 犏 ? xi yi 犏 犏 x2 y ? i i 犏 臌

法方程组为

?c1 n ? c2 ? x i ? c3 ? x i2 ? ? yi ? ? 2 3 ?c1 ? x i ? c2 ? x i ? c3 ? x i ? ? x i yi ? 2 3 4 2 c x ? c x ? c x ? x yi ? ? ? ? ? 1 i 2 i 3 i i ?

343

5.7.2 代数多项式拟合
3、一般情形
? ( x) ? c1?1 ( x) ? c2? 2 ( x) ? ?? cm? m ( x),(m ? n)

j 1 ( x) = 1, j 2 ( x) = x , j 3 ( x) = x2 , L , j m ( x) = xm- 1

超定方程组的 系数矩阵为

轾 1 犏 犏 1 犏 A= 犏 犏 犏 1 犏 臌
xi2 邋xi3 邋xi4 L 邋xim+ 1

x1 x2 L xn
L L L L L

x12 L
2 x2 L L L 2 xn L
m- 1 x 邋i xim xim+ 1

x1m- 1
m- 1 x2

,

( n > m)

m- 1 xn

轾 n 犏 犏 x 邋i 犏 AT A = 犏 邋xi2 犏 犏 犏 犏 m- 1 犏 邋xi 臌

邋xi xi2 xi3 L xim

,

xi2 m- 2

轾 yi 犏 犏 xy i i 犏 2 AT y = 犏 犏 xi yi 犏 M 犏 犏 m- 1 犏 臌 xi yi
344

5.7.2 代数多项式拟合
例 已知实验数据表,试用最小二乘法求经验公 式拟合这组数据。 解 作散点图,容易看出数据点接近一条直线, 因此设经验公式为 x 2 4 6 8 y = c0 + c1x, n=4, ? xi = 20 y 2 11 28 40 2

?

xi = 4 + 16 + 36 + 64 = 120

?

yi = 81

?

xi yi = 4 + 44 + 168 + 320 = 536

其法方程组为

ì 4c0 + 20c1 = 81 ? ? í ? ? ? 20c0 + 120c1 = 536

345

5.7.2 代数多项式拟合
解得

c0 = - 12.5, c1 = 6.55

得经验公式为 y= -12.5+6.55x
y 40 30 20 10 0
? ( 4 , 11)

(8, 40) ?

( 6, 28) ?

( 2, 2) ? 2 4 6 8 x
346

5.7.2 代数多项式拟合
例 已知一组观测数据表,试用最小二乘法求一 个多项式拟合这组数据。 x 0 1 2 3 4 5 y 5 2 1 1 2 3

解 作散点图,可以看出这些点接近一条抛物 线,因此设所求的多项式为

j ( x) = c0 + c1x + c2 x

2

347

5.7.2 代数多项式拟合
其法方程组为 ì 6c0 + 15c1 + 55c2 = 14 ? ? ? ? í 15c0 + 55c1 + 225c2 = 30 ? ? ? ? ? 55c0 + 225c1 + 979c2 = 122 得 c0=4.714 3 , c1= -2.785 7 , c2=0.500 0

j ( x) = 4.714 3 - 2.785 7 x + 0.500 0 x2
348

5.7.2 代数多项式拟合
x 0 1 y 5 2 2 1 3 1 4 2 5 3

349

5.7.2 代数多项式拟合
例 求一个经验函数 j ( x) = ae ,使它与观测数据拟合。(a, b 为常数 )
bx

x

1

2

3

4

5

6

7

8

y 14.3 20.5 27.4 36.6 49.1 64.6 87.8 117.6 解 对经验公式两边取对数得 ln j ( x) = ln a + bx 令 u=ln y , A=ln a , B=b
350

u=ln=A+Bx

5.7.2 代数多项式拟合
则函数 u=ln y(x) 的经验公式形为 u=ln=A+Bx 这是线性拟合,法方程组形为
ì ? ? An + B邋xi = í 2 ? A x + B x ? i = ? ? 邋i ui xiui
xi 2 = 204

n = 8,

?

xi = 36 ,

?

邋u =
i

ln yi = 29.978 7
351

?

xiui = 147.194 8 ,

5.7.2 代数多项式拟合
法方程组为
ì 8 A + 36B = 29.978 7 ? ? í ? ? ? 36 A + 204 B = 147.194 8

求得 A=2.430 482, B=0.292 6 ,

a = e = 11.36
j ( x) = 11.36e
0.292 6 x

A

352

5.7.2 代数多项式拟合

353


更多相关文档:

第二章 一元非线性方程的数值解法

第二章 一元非线性方程的数值解法 在科学和工程计算中,如电路和电力系统计算、 在科学和工程计算中,如电路和电力系统计算、非 线性微分和积分方程、非线性规划、...

实验报告二 一元非线性方程的解法..

实验报告二 一元非线性方程的解法.._数学_自然科学_专业资料。浙江大学城市学院...第二章 非线性方程(组)... 17页 免费 非线性方程的解法2-2 7页 免费 ...

计算方法 第2章 非线性方程数值解法

0 (2.1) 的数值解法,我们最为熟悉的非线性方程一元二次方程 ax2 ?计算方法讲稿 第二章 非线性方程数值解法 第二章本章将讨论非线性方程 非线性方程数值...

第二章 非线性方程(组)的数值解法

第三篇 第二章 非线性方程(组)的数值解法的 MATLAB 程序 第二章 2.1 2....2.9.2 求解 n 元非线性方程组的牛顿法及其 MATLAB 程序 解 n 元非线性...

6.6 一元非线性方程的解法习题课

第六章 一元非线性方程的解法习题课一、教学目标及基本要求通过对本节课的学习,使学生掌握方程求根的数值解法。 二、教学内容及学时分配本章主要介绍方程求根的...

一元非线性方程的数值解法

一元非线性方程的数值解法_理学_高等教育_教育专区。1 通过实习进一步掌握牛顿迭代...非线性方程(组)数值解法 34页 1下载券 第二章 非线性方程的数值... 21页 ...

第二章 非线性方程求根

第二章 非线性方程求根_数学_自然科学_专业资料。01...次代数方程以及某些特殊的高次方程或超越方程的解法...的一个不动点,也是方程(2.1)的一个根,如果 ? ...

实验四 一元非线性方程的迭代解法

实验四 一、 实验目的 一元非线性方程的迭代解法 通过该综合性实验的基本训练, 使学生进一步掌握 科学计算方法的基本知识; 掌握使用常用计算机语言编 程进行数值...

第二章 非线性方程(组)的数值解法

41页 免费 第18讲(非线性方程解法) 30页 1下载券第​二​章​ ​ ...[a,b] 的图形的 MATLAB 程序一 x=a:h:b; % h是步长 y=f(x); plot...

非线性方程的解法2-1

非线性方程的解法2-1_理学_高等教育_教育专区。第二章 非线性方程( 非线性方程...2第二章 非线性方程的数... 67页 1下载券 第二章 一元非线性方程... ...
更多相关标签:
非线性方程的数值解法 | 非线性方程数值解法 | 非线性方程的迭代解法 | 非线性方程解法 | 非线性方程的解法 | 多元非线性方程 | 一元非线性方程 | 多元非线性方程求解 |
网站地图

文档资料共享网 nexoncn.com copyright ©right 2010-2020。
文档资料共享网内容来自网络,如有侵犯请联系客服。email:zhit325@126.com