Transcript 第三章
第3章
函数逼近与快速傅里叶变换
1
3.1
函数逼近的基本概念
3.2
正交多项式
3.3
最佳平方逼近
3.4
曲线拟合的最小二乘法
3.5
有理逼近
3.6
三角多项式与快速傅里叶变换
3.1
3.1.1
问题
函数逼近的基本概念
函数逼近与函数空间
1、数值计算中经常要计算函数值,如计算机中计算
基本初等函数及其他特殊函数;
2、当函数只在有限点集上给定函数值,要在包含该
点集的区间上用公式给出函数的简单表达式.
这些都涉及到在区间 [a, b] 上用简单函数逼近已知复杂
函数的问题,这就是函数逼近问题.
2
插值法就是函数逼近问题的一种.
本章讨论的函数逼近,是指“对函数类 A 中给定的函数
f (x), 记作 f ( x) A , 要在另一类简单的便于计算的函数类
B 中求函数 p( x) B ,使 p (x) 与 f (x) 的误差在某种度量
意义下最小”.
函数类 A通常是区间 [a, b] 上的连续函数,记作 C[a, b] ,
称为连续函数空间.
3
函数类 B通常为 n 次多项式,有理函数或分段低次多项
式等.
数学上常把在各种集合中引入某些不同的确定关系称为
赋予集合以某种空间结构,并将这样的集合称为空间.
例如将所有实 n维向量组成的集合,按向量加法及向量
与数的乘法构成实数域上的线性空间, 记作 R n, 称为 n 维
向量空间.
4
对次数不超过 n( n为正整数)的实系数多项式全体,
按通常多项式与多项式加法及数与多项式乘法也构成数域
R上一个线性空间,用 H n表示,称为多项式空间.
所有定义在 [a, b] 上的连续函数集合,按函数加法和
数与函数乘法构成数域 R 上的线性空间,记作 C[a, b].
类似地, 记 C p [a, b]为具有 p 阶连续导数的函数空间.
5
定义1
设集合 S 是数域 P上的线性空间,元素
x1 ,, xn S , 如果存在不全为零的数 1 ,, n P ,
使得
1 x1 n xn 0,
(1.1)
则称 x1 ,, xn 线性相关.
否则,若等式(1.1)只对 1 2 n 0成立,
则称 x1 , , xn线性无关.
6
若线性空间 S 是由 n 个线性无关元素 x1 ,, xn 生成的,
x S 都有
x 1 x1 n xn
则 x1 ,, xn 称为空间 S 的一组基,记为
S span{x1 ,, xn }
并称空间 S 为 n 维空间,系数 1 ,, n 称为 x 在基
x1 ,, xn下的坐标, 记作 (1 ,, n ).
如果 S 中有无限个线性无关元素 x1 , , xn ,, 则称 S
为无限维线性空间.
7
考察次数不超过 n次的多项式集合 H n, 其元素
p( x) H n 表示为
p( x) a0 a1 x an x n ,
(1.2)
它由 n 1个系数 (a0 , a1 ,, an ) 唯一确定.
1, x, , x n是线性无关的,它是 H n 的一组基,故
H n span{1, x,, x n },
且 (a0 , a1 ,, an ) 是 p (x) 的坐标向量,H n 是 n 1维的.
8
对连续函数 f ( x) C[a, b],它不能用有限个线性无关的
函数表示,故 C[a, b] 是无限维的,但它的任一元素 f (x)
均可用有限维的 p( x) H n 逼近,使误差
max f ( x) p( x)
a x b
( 为任给的小正数),这就是著名的魏尔斯特拉斯定理.
9
定理1
设 f ( x) C[a, b] , 则对任何 0 ,总存在一
个代数多项式 p (x) , 使
f ( x) p ( x )
在 [a, b] 上一致成立.
伯恩斯坦1912年给出的证明是一种构造性证明.
他根据函数整体逼近的特性构造出伯恩斯坦多项式
n
k
Bn ( f , x) f ( ) Pk ( x)
n
k 0
(1.3)
10
其中
n k
Pk ( x) x (1 x) n k ,
k
n n(n 1) (n k 1)
k!
k
为二项式展开系数,并证明了
lim Bn ( f , x) f ( x) 在 [0,1] 上一致成立;
n
若 f (x) 在 [0,1]上 m阶导数连续,则
lim Bn( m ) ( f , x) f
n
( m)
( x).
这个结果不但证明了定理1,而且由(1.3)给出了 f (x)
的一个逼近多项式,但它收敛太慢,实际中很少使用.
11
更一般地,可用一组在 C[a, b] 上线性无关的函数集合
i ( x) in0 来逼近 f ( x) C[a, b],此时元素
( x) span{0 ( x), 1 ( x), , n ( x)} C[a, b]
可表示为
( x) a00 ( x) a11 ( x) ann ( x).
(1.4)
函数逼近问题就是对任何 f ( x) C[a, b],在子空间Φ
*
*
找一个元素 ( x) , 使 f ( x) ( x)在某种意义下最小.
12
3.1.2
范数与赋范线性空间
为了对线性空间中元素大小进行衡量,需要引进范数
定义,它是 R n空间中向量长度概念的直接推广.
13
定义2
设 S为线性空间,x S ,若存在唯一实数‖·‖,
满足条件:
(1)
x 0, 当且仅当 x 0 时, x 0; (正定性)
(2)
x x ,
(3)
x y x y ,
R;
(齐次性)
x, y S . (三角不等式)
则称‖·‖为线性空间 S上的范数,S 与‖·‖一起称为赋范
线性空间,记为 X .
14
例如,在 R n上的向量 x ( x1 , , xn )T R n , 三种常
用范数为
x
x
1
max xi ,
1i n
n
i 1
2
称为 1-范数,
xi ,
n
x
称为 范数或最大范数,
1
2
( xi2 ) ,
称为 2-范数.
i 1
15
实际上任何向量的实值函数,只要满足上述三个条件,
就可以定义成一种向量范数.
在 R 中,满足‖·‖2 =1 ,即 x1 x2
2
2
2
1的向量
为单位圆,
满足‖·‖∞ =1 ,即 max{ x1 , x2 } 1 的向量为单位正
方形,
而满足‖·‖1 =1 的向量 x1 x2 1 则为对角线长
度为1的菱形.
16
所以说,范数是对向量长度的度量,度量方式不同,
结果也不一样,但不同范数之间是存在等价关系的.
17
类似地,对连续函数空间 C[a, b] ,若 f ( x) C[a, b] ,
f
f
f
1
a x b
b
称为 1-范数,
f ( x) dx,
a
b
2
称为 范数,
max f ( x) ,
1
2
( f 2 ( x) dx) ,
称为 2-范数.
a
可以验证这样定义的范数均满足定义2中的三个条件.
18
3.1.3
内积与内积空间
在线性代数中, R n中两个向量 x ( x1 , , xn )T 及
y ( y1 , , yn )T 的内积定义为
( x, y) x1 y1 xn yn .
(1.5)
若将它推广到一般的线性空间 X,则有下面的定义.
19
定义3 X 是数域K(R或C)上的线性空间,对
u, v X, 有K中一个数与之对应,记为 (u, v) ,它满足
以下条件:
(1)
(u, v) (v, u ),
u, v X;
(2)
(u, v) (u, v),
(3)
(u v, w) (u, w) (v, w),
(4)
(u, u) 0,
K, u, v X;
u, v, w X;
当且仅当u 0 时,
(u, u) 0.
则称 (u, v)为X上 u 与 v 的内积.
20
定义了内积的线性空间称为内积空间.
定义中(1)的右端 (u , v) 称为 (u, v) 的共轭,
当K为实数域R时 (u, v) (v, u ).
如果 (u, v) 0 ,则称 u与 v 正交,这是向量相互垂
直概念的推广.
21
定理2 设X为一个内积空间, 对 u , v X , 有
(u, v) (u, u )(v, v).
2
(1.6)
称为柯西-施瓦茨(Cauchy-Schwarz)不等式.
证明
当 v 0 时(1.6)式显然成立.
现设 v 0 ,则 (v, v) 0, 且对任何数 有
0 (u v, u v) (u, u ) 2 (u, v) 2 (v, v).
取 (u, v) /( v, v) ,代入上式右端,得
22
(u, u ) 2
(u , v)
2
(v, v)
(u, v)
2
(v, v)
0
即得 v 0时
(u, v) (u, u )(v, v).
2
23
定理3
设X为一个内积空间,u1 , u2 , , un X,
矩阵
(u1 , u1 )
(u , u )
G 1 2
(u1 , un )
(u2 , u1 )
(u2 , u2 )
(u 2 , un )
(u n , u1 )
(u n , u2 )
(u n , u n )
(1.7)
称为格拉姆(Gram)
则 G 非奇异的充分必要条件是 u1 , u2 ,, un 线性无关.
24
G非奇异等价于 det G 0,其充要条件是齐次
证明
方程组
n
n
j 1
j 1
( j u j , uk ) (u j , uk ) j 0, k 1,2,, n(1.8)
只有零解; 而
n
j 1
j
u j 1u1 2u2 nun 0
n
n
j 1
j 0
(1.9)
( j u j , j u j ) 0
n
( j u j , uk ) 0,
j 0
k 1,2,, n
25
从以上等价关系知,det G 0 等价于从(1.8)推出
1 2 n 0,
n
n
而后者等价于从(1.9)推出
n
0,
( j u j , uk ) (u j , uk ) j
1 0, k2
(1.8)
1,2,, n
j 1
j 1
即 u1 , u2 ,, un 线性无关.
n
在内积空间X
u u
j 1
j
j
1 1
2u 2 n u n 0
u X, 记
u
(u, u )
(1.10)
26
利用
(u v) u
2
2
2 u v v
2
(u, u ) 2(u, v) (v, v)
(u v, u v) u v
2
两端开方即得三角不等式
uv u v
(1.11)
27
例1
设
R n 与 C n的内积.
T
x, y R n , x ( x1 , , xn )T , y ( y1 ,, yn ) ,
n
( x, y ) xi yi
(1.12)
i 1
向量2-范数为
x
1
2
2
n
( x, x) ( xi2 )
1
2
i 1
28
若给定实数 i 0(i 1,2,, n), 称{i }为权系数,
R n 上的加权内积为
n
( x, y ) i xi yi
(1.13)
i 1
相应的范数为
n
x
2
1
2
( i xi2 ) .
i 1
当 i 1(i 1,2,, n) 时,
(1.13)就是前面定义的内积.
29
如果 x, y C n,带权内积定义为
n
( x, y ) i xi yi
(1.14)
i 1
这里 {i } 仍为正实数序列,yi 为 yi 的共轭.
在 C[a, b]上也可以类似定义带权内积.
30
定义4
设 [a, b] 是有限或无限区间,在 [a, b] 上的非负
函数 (x)满足条件:
(1) a x k ( x)dx 存在且为有限值 (k 0,1, ),
b
(2) 对 [a, b] 上的非负连续函数 g (x) ,如果
b
a
g ( x) ( x)dx 0,
g ( x) 0.
则称 (x) 为[a, b]上的一个权函数.
31
例2 C[a, b] 上的内积.
设 f ( x), g ( x) C[a, b], (x) 是 [a, b] 上给定的权函数,
则可定义内积
b
( f ( x), g ( x)) ( x) f ( x) g ( x)dx.
(1.15)
a
由此内积导出的范数为
f ( x)
1
2
2
1
2
b
( f ( x), f ( x)) ( x) f 2 ( x)dx (1.16)
a
称(1.15)和(1.16)为带权 (x) 的内积和范数.
32
常用的是 ( x) 1 的情形,即
b
( f ( x), g ( x)) f ( x) g ( x)dx.
a
f ( x)
2
b
f 2 ( x)dx
a
1
2
33
若 0 , 1 ,, n是 C[a, b]中的线性无关函数族,记
span{0 , 1 ,, n }, 它的格拉姆矩阵为
G G(0 , 1 , , n )
( 0 , 0 )
( , )
1 0
( n , 0 )
( 0 , 1 )
(1 , 1 )
( n , 1 )
( 0 , n )
(1 , n )
( n , n )
(1.17)
根据定理3可知 0 , 1 ,, n 线性无关的充要条件是
det G(0 , 1 , , n ) 0.
34
3.1.4
最佳逼近
函数逼近主要讨论给定 f C[a, b],求它的最佳逼近
多项式的问题.
若 P* ( x) H n使误差
f ( x) P* ( x) min f ( x) P( x)
PH n
则称 P* ( x) 是 f (x)在 [ a, b]上的最佳逼近多项式.
若 P( x) span{0 , 1 ,, n }则称相应的 P* ( x)
为最佳逼近函数.
通常将范数 取为
或 2.
35
若取 ,即
f ( x) P * ( x)
min f ( x) P( x)
PH n
(1.18)
min max f ( x) P( x) ,
PH n a x b
则称 P* ( x) 是 f (x)在 [ a, b]上的最优一致逼近多项式.
求 P* ( x)就是求 [ a, b]上使最大误差 max f ( x) P( x)
a x b
最小的多项式.
36
若取 ,即
2
f ( x) P * ( x)
2
2
min f ( x) P( x)
PH n
min
PH n
b
a
2
2
(1.19)
[ f ( x) P( x)]2 dx,
则称 P* ( x) 是 f (x)在 [ a, b]上的最佳平方逼近多项式.
若 f (x)是 [ a, b]上的一个列表函数,在
a x0 x1 xm b
上给出 f ( xi )(i 0,1,, m) ,要求 P* 使
f P * 2 min f P
P
2
min
P
m
2
[
f
(
x
)
P
(
x
)]
i
i
(1.20)
i 0
则称 P* ( x) 为 f (x)的最小二乘拟合.
37
3.2
正交多项式
3.2.1
定义5
正交函数族与正交多项式
若 f ( x), g ( x) C[a, b], ( x )为 [ a, b]
上的权函数且满足
( f ( x), g ( x))
则称 f
b
a
( x) f ( x) g ( x)dx 0.
(2.1)
( x )与 g ( x )在 [ a, b]上带权 ( x )正交.
38
若函数族 0 ( x), 1 ( x), , n ( x), 满足关系
( j , k )
b
a
j k.
0,
( x) j ( x) k ( x)dx
j k.
Ak 0,
(2.2)
则称{ k ( x)}是 [a, b] 上带权 (x)的正交函数族.
若 Ak 1,则称之为标准正交函数族.
39
三角函数族
1, cos x, sin x, cos 2 x, sin 2 x,
就是在区间 [ , ] 上的正交函数族.
定义6
设 n ( x)是 [a, b] 上首项系数 an 0 的 n次多
项式, (x) 为 [ a, b] 上权函数,如果多项式序列 { n ( x)}
0
满足关系式(2.2),则称多项式序列 { n ( x)}0 为在 [a, b] 上
带权 (x) 正交,称 n ( x) 为 [a, b] 上带权 (x) 的 n次正
交多项式.
40
只要给定区间 [a, b] 及权函数 (x) ,均可由一族线性
无关的幂函数 {1, x,, x n ,}, 利用逐个正交化手续构造
出正交多项式序列 { n ( x)}0 :
0 ( x) 1,
n 1
n ( x) x n
j 0
( x n , j ( x))
j ( x)
( j ( x), j ( x))
(n 1,2, ).
(2.3)
41
正交多项式 n ( x)的最高次项系数为1. 而若 { n ( x)}0
是正交多项式,则 0 ( x), 1 ( x), , n ( x) 在 [ a, b]上是
线性无关的.
事实上,若
c00 ( x) c11 ( x) cn n ( x) 0,
用 ( x) j ( x)( j 0,1, , n)乘上式并积分得
b
b
c0 ( x ) 0 ( x ) j ( x )dx c1 ( x )1 ( x ) j ( x )dx
a
a
b
c j ( x ) j ( x ) j ( x )dx
a
b
cn ( x ) n ( x ) j ( x )dx 0,
a
42
利用正交性有
b
c j ( x) j ( x) j ( x)dx 0.
a
b
由于 a ( x) j ( x) j ( x)dx 0,故 c j 0( j 0,1, , n)
即0 ( x), 1 ( x), , n ( x) 线性无关.
关于正交多项式,有
(1) 任何 P( x) H n 均可表示为 0 ( x), 1 ( x), , n ( x)
的线性组合. 即
n
P( x) c j j ( x).
j 0
43
(2) n ( x)与任一次数小于 n的多项式 P( x) H n 1
正交. 即
( n ( x), P)
b
a
( x) n ( x) P( x)dx 0.
除此之外,还有
44
定理4
设 { n ( x)}0 是 [ a, b]上带权 ( x )的正交多项式,
对 n 0成立关系
n1 ( x) ( x n )n ( x) nn1 ( x)
(n 0,1, ).
(2.4)
其中
0 ( x) 1,
1 ( x) 0,
n ( x n ( x), n ( x)) /( n ( x), n ( x)),
n ( n ( x), n ( x)) /( n 1 ( x), n 1 ( x)),
(n 1,2, ).
b
2
这里 ( x n ( x), n ( x)) a xn ( x) ( x)dx.
45
设 { n ( x)}0 是 [ a, b]上带权 ( x )的正交多项式,
则 n ( x)( n 1)在区间 (a, b)内有n 个不同的零点.
定理5
证明
假定 n (x) 在(a, b) 内的零点都是偶数重的,则
n (x)在 [ a, b]符号保持不变,这与
( n ( x), 0 ( x))
b
a
( x) n ( x)0 ( x)dx 0.
矛盾.故 n (x) 在 (a, b)内的零点不可能全是偶数重的,现设
xi (i 1,2,, l )为 n (x)在 (a, b)内的奇数重零点,不妨设
a x1 x2 xl b,
则 n (x) 在 xi (i 1,2,, l ) 处变号.
46
令
q( x) ( x x1 )( x x2 ) ( x xl ),
于是 n ( x)q( x)在 [ a, b]上不变号,则得
( n , q)
b
a
( x) n ( x)q( x)dx 0.
若 l n,由 { n ( x)}0 的正交性可知
( n , q)
b
a
( x) n ( x)q( x)dx 0,
这与 ( n , q) 0矛盾,故 l n.而 n (x)只有 n个零点,
故 l n ,即 n个零点都是单重的.
47
3.2.2
勒让德多项式
当区间为[1, 1] ,权函数 ( x) 1时,由 {1, x,, x n ,}
正交化得到的多项式就称为勒让德(Legendre)多项式,
并用 P0 ( x), P1 ( x), , Pn ( x), 表示.
罗德利克(Rodrigul)给出了简单的表达式
P0 ( x) 1,
1 dn
2
n
Pn ( x) n
{(
x
1
)
}
n
2 n! dx
(n 1,2, ),
(2.5)
48
由于 ( x 2 1) n是 2n 次多项式,所以对其求 n 阶导数后得
Pn ( x)
1
n
n 1
(
2
n
)(
2
n
1
)
(
n
1
)
x
a
x
a0 ,
n 1
n
2 n!
于是得首项 x 的系数 an
n
( 2n)!
.
n
2
2 ( n!)
最高项系数为1的勒让德多项式为
~
Pn ( x)
n! d n
2
n
[(
x
1
)
].
n
(2n)! dx
(2.6)
49
勒让德多项式重要性质:
性质1
正交性
0,
1 Pn ( x) Pm ( x)dx 2 ,
2n 1
1
证明
m n;
m n.
(2.7)
令 ( x) ( x 2 1) n,则
( k ) (1) 0
(k 0,1, , n 1).
设 Q(x)是在区间 [1, 1] 上 n 阶连续可微的函数,由分部
积分知
50
1
1
Pn ( x)Q( x)dx
1
1
(n)
Q
(
x
)
( x)dx
n
1
2 n!
1
1
n Q( x) ( n 1) ( x)dx
2 n! 1
( 1) n
n
2 n!
1
1
Q ( n ) ( x ) ( x ) dx
下面分两种情况讨论:
(1) 若 Q(x) 是次数小于 n的多项式,则 Q ( n ) ( x) 0,
故得
1
1
Pn ( x) Pm ( x)dx 0, 当n m.
51
(2) 若
则
于是
Q( x) Pn ( x)
(2n)! n
1
(n)
x ,
(
x
)
n
2
n
2 (n!)
2 n!
Q ( n ) ( x) Pn( n ) ( x)
( 2n)!
,
n
2 n!
n
(
1
)
(2n)! 1
2
2
n
P
(
x
)
dx
(
x
1
)
dx
n
2
2
1
1
2 (n!)
1
(2n)!
2n
2 (n!) 2
1
1
(1 x 2 ) n dx.
由于
1
0
故
(1 x ) dx cos 2 n 1 tdt
2
n
2
0
1
1
Pn2 ( x) dx
2
.
2n 1
52
性质2
奇偶性
Pn ( x) (1) n Pn ( x).
(2.8)
由于 ( x) ( x 2 1) n 是偶次多项式,经过偶次求导仍为
偶次多项式,经过奇次求导则为奇次多项式,故 n为偶数时
Pn ( x )为偶函数,n为奇数时 Pn (x)为奇函数,于是(2.8)成
立.
53
性质3
递推关系
考虑 n 1次多项式 xPn ( x), 它可表示为
xPn ( x) a0 P0 ( x) a1 P1 ( x) an 1 Pn 1 ( x).
两边乘 Pk ( x), 并从-1到1积分,得
1
1
1
xPn ( x) Pk ( x)dx ak Pk2 ( x)dx.
1
当 k n 2 时, xPk ( x)次数小于等于 n 1,上式左端积分
为0, 故得 ak 0.
54
当 k n 时,xPn2 ( x)为奇函数,左端积分仍为0, 故
an 0.
于是
其中
xPn ( x) an 1 Pn 1 ( x) an 1 Pn 1 ( x)
an 1
an 1
2n 1 1
xPn ( x ) Pn 1 ( x ) dx
1
2
2n 1
2n
n
,
2
2
4n 1
2n 1
2n 3 1
xPn ( x ) Pn 1 ( x) dx
1
2
2n 3
2( n 1)
n 1
,
2
( 2n 1)( 2n 3)
2n 1
55
从而得到以下的递推公式
(n 1) Pn 1 ( x) (2n 1) xPn ( x) nPn 1 ( x)
(2.9)
(n 1,2, ),
由 P0 ( x) 1, P1 ( x) x, 利用上述递推公式就可推出
P2 ( x) (3x 2 1) / 2,
P3 ( x) (5 x 3 3 x) / 2,
P4 ( x) (35 x 4 30 x 2 3) / 8,
P5 ( x) (63 x 5 70 x 3 15 x) / 8,
P6 ( x) (231x 6 315 x 4 105 x 2 5) / 16,
56
图3-1给出了P0 ( x), P1 ( x), P2 ( x), P3 ( x) 的图形.
图3-1
57
性质4
Pn ( x )在区间 [1, 1] 内有 n 个不同的实零点.
58
3.2.3
切比雪夫多项式
当权函数 ( x)
1
1 x2
,区间为 [1, 1]时,由序
列 {1, x,, x n ,} 正交化得到的正交多项式就是切比雪夫
(Chebyshev)多项式.
它可表示为
Tn ( x) cos( n arccos x),
x 1.
(2.10)
若令 x cos , 则 Tn ( x) cos n , 0 .
59
切比雪夫多项式有很多重要性质:
性质1
递推关系
Tn 1 ( x) 2 xTn ( x) Tn 1 ( x)
T0 ( x) 1,
(n 1,2, ),
T1 ( x) x.
(2.11)
这只要在三角恒等式
cos( n 1) 2 cos cos n cos( n 1)
(n 1)
中,令 x cos 即得.
60
由(2.11)可推出
T2 ( x) 2 x 2 1,
3
(
x
)
4
x
3 x,
TT
(
x
)
2
xT
3
n 1
n ( x) Tn 1 ( x)
T4 ( x) 8 x 4 8 x 2 1,
T5 ( x) 16 x 5 20 x 3 5 x,
T6 ( x) 32 x 6 48 x 4 18 x 2 1,
T0 ( x), T1 ( x), T2 ( x), T3 ( x) 的函数图形见图3-2.
61
图3-2
由递推关系(2.11)还可得到 Tn (x) 的最高次项系数是
2 n 1.
Tn 1 ( x) 2 xTn ( x) Tn 1 ( x)
62
性质2
切比雪夫多项式 {Tk ( x)} 在区间 [1, 1]上带权
( x) 1 / 1 x 2 正交,且
0,
1 T ( x )T ( x ) dx
n
m
,
1
2
1 x
2
,
n m;
n m 0;
(2.12)
n m 0.
令 x cos ,则 dx sin d ,于是
0, n m;
1 T ( x )T ( x ) dx
n
m
cos n cos md
, n m 0;
1
2
0
1 x
2
, n m 0.
63
性质3
T2 k ( x) 只含 x 的偶次幂,T2 k 1 ( x)只含 x 的奇次幂.
这个性质由递推关系直接得到.
性质4
Tn (x) 在区间 [1, 1]上有 n 个零点
2k 1
xk cos
,
2n
性质5
k 1,2,, n.
Tn (x) 的首项 x n 的系数为 2 n 1.
~
若令 T~0 ( x) 1, T~n ( x) 1
,则
T
(
x
),
n
1
,
2
,
Tn ( x)
n 1 n
2
是首项系数为1的切比雪夫多项式.
64
~ 为所有次数小于等于 的首项系数为1的多项
若记 H
n
n
式集合,对 T~n ( x)有以下性质:
定理6
设 T~n ( x)是首项系数为1的切比雪夫多项式,则
~
max Tn ( x) max P( x) ,
1 x 1
1 x 1
~
P( x) H n .
且
~
max Tn ( x)
1 x 1
1
.
n 1
2
65
~ 中,
定理6表明在所有首项系数为1的 n次多项式集合 H
n
~
~
~ 中最大值最小的多项
所以
是
Tn ( x) min
P
(
x
)
,
T
(
x
)
H
~
n
n
PH
n
式,即
~
max Tn ( x) min
~ max P ( x )
1 x 1
PH n 1 x 1
1
.
n 1
2
(2.13)
利用这一结论,可求 P( x) H n 在 H n 1中的最佳(一致)
逼近多项式.
66
例3 求 f ( x) 2 x3 x 2 2 x 1 在 [1, 1]上的最佳2次逼
近多项式.
解 由题意,所求最佳逼近多项式 P2* ( x) 应满足
max f ( x) P2* ( x) min .
1 x 1
由定理6可知, 当
1
3
3
f ( x) P ( x) T3 ( x) 2 x x
2
2
*
2
*
时,多项式 f ( x) P2 ( x)
故
67
1
7
2
P ( x) f ( x) T3 ( x) x x 1
2
2
*
2
就是 f (x) 在 [1, 1]上的最佳2次逼近多项式.
由于切比雪夫多项式是在区间 [1, 1]上定义的,对于
一般区间 [ a, b],要通过变量替换变换到 [1, 1] ,可令
1
x [(b a) t a b],
2
则可将 x [a, b] 变换到 t [1, 1].
(2.14)
68
3.2.4
切比雪夫多项式零点插值
切比雪夫多项式 Tn (x)在区间 [1, 1]上有 n 个零点
xk cos
2k 1
,
2n
k 1,2,, n.
和 n 1个极值点(包括端点)
k
xk cos
,
n
k 0,1, , n.
这两组点称为切比雪夫点,它们在插值中有重要作用.
69
从图3-3可以看到切比雪夫点恰好是单位圆周上等距分
布点的横坐标,这些点的横坐标在接近区间 [1, 1]的端点
处是密集的.
利用切比雪夫点做插值,可
使插值区间最大误差最小化.
设插值点 x0 , x1 ,, xn [1,1]
f C n 1[1,1], Ln ( x)为相应的 n
次拉格朗日插值多项式,则余项
图3-3
f ( n1 ( )
Rn ( x) f ( x) Ln ( x)
n1 ( x),
(n 1)!
70
于是
max f ( x) Ln ( x)
1 x 1
M n 1
max ( x x0 )( x x1 ) ( x xn ) ,
(n 1)! 1 x 1
其中
M n 1 f ( n 1) ( x)
max f ( n 1) ( x)
1 x 1
是由被插函数决定的.
如果插值节点为 Tn1 ( x)的零点
xk cos
2k 1
,
2(n 1)
k 0,1, , n,
71
则由(2.13)可得
1
~
max n 1 ( x) max Tn 1 ( x) n .
1 x 1
1 x 1
2
由此可导出插值误差最小化的理论.
定理7
设插值节点 x0 , x1 ,, xn 为切比雪夫多项式 Tn1 ( x)
的零点,被插函数 f C n 1[1,1], Ln ( x) 为相应的插值多项
式,则
1
max f ( x) Ln ( x) n
f ( n 1) ( x) .
1 x 1
2 (n 1)!
(2.15)
72
对于一般区间 [ a, b]上的插值只要利用变换(2.14)则
可得到相应结果,此时插值节点为
ba
2k 1
ab
xk
cos
,
2
2(n 1)
2
k 0,1, , n,
73
例4
求 f ( x) e x在 [0, 1]上的四次拉格朗日插值多项式
x
max
e
L4 ( x) .
,插值节点用
的零点并估计误差
L4 ( x)
T5 ( x)
0 x 1
解
利用 T5 ( x)的零点和区间变换可知节点
xk
1
2k 1
(1 cos
),
2
10
k 0,1,2,3,4,
即
x0 0.97553,
x1 0.79390,
x3 0.20611,
x4 0.02447,
x2 0.5,
对应的拉格朗日插值多项式为
74
L4 ( x) 1.000 022 74 0.998 862 33x 0.509 022 51x 2
0.141 841 05 x 3 0.068 494 35 x 4
利用(2.15)可得误差估计
M n 1 (b a) n 1
max e L4 ( x)
,
2 n 1
0 x 1
(n 1)! 2
x
n 4,
而
M n1 f (5) ( x)
ex
e1 2.72,
于是有
max e x L4 ( x)
0 x 1
e 1
2.72 1
. 9
4.4 10 5.
5! 2
6 10240
75
在第2章中曾经指出,高次插值会出现龙格现象,一般
Ln (x)不收敛于 f (x),因此并不适用. 但若用切比雪夫多项
式零点插值却可避免龙格现象,可保证整个区间上收敛.
76
1
,在 [5,5]上利用T11 ( x)的零点
2
1 x
~
作插值点,构造10次拉格朗日插值多项式 L
. 与第2章
10 ( x)
例5
设 f ( x)
得到的等距节点造出的 L10 ( x)近似 f (x) 作比较.
解
在[1, 1] 上的10次切比雪夫多项式 T11 ( x) 的零点为
t k cos
21 2k
,
22
k 0,1, ,10.
作变换 xk 5tk ,
k 0,1,,10.它们是 (5,5)内的插值点,
~
由此得到 y f (x) 在 [5,5]上的拉格朗日插值多项式 L
10 ( x ).
77
~
~
的图形见图3-4,从图上可见
L10 ( x)
f ( x), L10 ( x), L10 ( x)
没有出现龙格现象.
图3-4
78
3.2.5 其他常用的正交多项式
区间 [ a, b]及权函数 (x)不同,则得到的正交多项式也
不同.
除上述两种最重要的正交多项式外,下面是三种较常
用的正交多项式.
1. 第二类切比雪夫多项式
在区间[1, 1] 上带权 ( x)
1 x 2 的正交多项式
称为第二类切比雪夫多项式.
79
表达式为
U n ( x)
sin[( n 1) arccos x]
1 x2
.
(2.16)
令 x cos ,可得
1
U n ( x)U m ( x) 1 x dx sin( n 1) sin( m 1)d
1
2
0
0, m n;
, m n,
2
即 {U n ( x)}是 [1, 1]上带权 ( x)
1 x 2 的正交多项式族.
80
递推关系
U 0 ( x) 1, U1 ( x) 2 x,
U n1 ( x) 2 xUn ( x) U n 1 ( x),
(n 1,2, ).
81
2. 拉盖尔多项式
在区间 [0,)上带权 ( x) e x 的正交多项式称为拉
盖尔(Laguerre)多项式.
其表达式为
n
d
n x
Ln ( x) e x
(
x
e ).
n
dx
(2.17)
正交性质
0
0, m n;
e Ln ( x) Lm ( x)dx
2
(n!) , m n,
x
82
递推关系
L0 ( x) 1,
L1 ( x) 1 x,
Ln 1 ( x) (1 2n x) Ln ( x) n 2 Ln 1 ( x),
(n 1,2, ).
83
3. 埃尔米特多项式
在区间 (,) 上带权 ( x) e
x2
的正交多项式称
为埃尔米特多项式.
表达式
H n ( x) (1) e
n
x2
dn
x2
(e ),
n
dx
(2.18)
正交关系
e
x2
0, m n;
H n ( x) H m ( x)dx n
2 n! , m n,
84
递推关系
H 0 ( x) 1,
H1 ( x) 2 x,
H n1 ( x) 2 xHn ( x) 2nH n 1 ( x),
(n 1,2, ).
85
3.3
最佳平方逼近
86
最佳平方逼近及其计算
3.3.1
对 f ( x) C[a, b] 及 C[a, b] 中的一个子集
span{0 ( x), 1 ( x), , n ( x)}
若存在 S * ( x) ,使
f ( x) S ( x)
*
2
2
min f ( x) S ( x)
S ( x )
min
S ( x )
b
a
2
(3.1)
2
( x)[ f ( x) S ( x)]2 dx.
则称 S * ( x) 是 f (x) 在子集 C[a, b]中的最佳平方逼近
函数.
87
由(3.1)可知该问题等价于求多元函数
n
I (a0 , a1 , , an ) ( x)[ a j j ( x) f ( x)]2 dx (3.2)
a
b
j 0
的最小值.
I (a0 , a1 , , an ) 是关于 a0 , a1 , , an 的二次函数,
利用多元函数求极值的必要条件
I
0
ak
即
(k 0,1,, n),
n
b
I
2 ( x)[ a j j ( x) f ( x)] k ( x)dx
a
ak
j 0
(k 0,1,, n),
88
于是有
n
(
j 0
k
( x), j ( x)) a j ( f ( x), k ( x)) (k 0,1, , n). (3.3)
这个关于 a0 , a1 ,, an 的线性方程组,称为法方程.
由于 0 ( x), 1 ( x), , n ( x) 线性无关,故
det G(0 , 1 ,, n ) 0
于是方程组(4.3)有唯一解 ak ak*
(k 0,1, , n),
S * ( x) a0*0 ( x) an* n ( x).
89
下面证明 S * ( x)满足(3.1),即对任何 S ( x) , 有
b
a
b
( x)[ f ( x) S ( x)] dx ( x)[ f ( x) S ( x)]2 dx
*
2
a
(3.4)
为此只要考虑
2
2
(3.1)
f ( x) S ( x) min2 f ( x) S ( x) 2
2 S (Sx
D ( x)[ f ( x)
)])dx ( x)[ f ( x) S * ( x)]2 dx
(x
*
b
a
b
a
b
b
a
*min 2 ( x)[ f ( x) S ( x)]2 dx.
( x)[ f ( x) S S(( xx ))]
a
dx
b
2 ( x)[ S * ( x) S ( x)][ f ( x) S * ( x)]dx.
a
90
由于S * ( x)的系数 a k* 是方程(3.3)的解,故
b
a
( x)[ f ( x) S * ( x)]k ( x)dx 0 (k 0,1,, n),
n
从而上式第二个积分为0,
x), k ( x)) (k 0,1, , n). (3.3)
(k ( x), j ( x)) a j ( f (于是
j 0
b
D ( x)[ S ( x) S * ( x)]2 dx 0,
a
故(3.4)成立.
这就证明了S * ( x)是 f (x) 在 中的最佳平方逼近函数.
b
b
*
2
2
(
x
)[
f
(
x
)
S
(
x
)]
dx
(
x
)[
f
(
x
)
S
(
x
)]
dx
a
a
(3.4)
91
若令 ( x) f ( x) S * ( x), 则平方误差为
( x)
2
2
( f ( x) S * ( x), f ( x) S * ( x))
( f ( x), f ( x)) ( S * ( x), f ( x))
f ( x)
2
2
n
ak* ( k ( x), f ( x)).
(3.5)
k 0
若取 k ( x) x k , ( x) 1, f ( x) C[0, 1], 则要在 H n
中求 n次最佳平方逼近多项式
S * ( x) a0* a1* x an* x n ,
92
此时
( j ( x), k ( x))
( f ( x), k ( x))
1
0
1
0
x k j dx
1
,
k j 1
f ( x) x k dx d k .
n
若用 H表示 Gn G (1, x, , x )
1
1/ 2
H
1 /( n 1)
1/ 2
1/ 3
1 /( n 2)
1 /( n 1)
1 /( n 2)
1 /( 2n 1)
即
(3.6)
称为希尔伯特(Hilbert)矩阵.
93
T
T
记 a (a0 , a1 ,, an ) , d (d 0 , d1 ,, d n ) , 则
Ha d
(3.7)
*
a
a
的解 k
k ( k 0,1, , n) 即为所求.
94
设 f ( x) 1 x 2 , 求 [0, 1]上的一次最佳平方
例6
逼近多项式.
解 利用(3.7),得
1
d0
0
1
1 x dx ln( 1
2
2
2
2)
1.147,
2
2 2 1
1
2 3/ 2
2
0.609,
d1 x 1 x dx (1 x )
0
3
3
0
1
1
得方程组
1 / 2 a0 1.147
1
1 / 2 1 / 3 a 0.609,
1
95
解之
a0 0.934,
a1 0.426,
故
S1* ( x) 0.934 0.426 x.
平方误差
( x)
2
2
( f ( x), f ( x)) ( S1* ( x), f ( x))
1
(1 x 2 )dx 0.426d1 0.934d 0 0.0026.
0
最大误差
( x)
max
0 x 1
1 x 2 S1* ( x) 0.066.
96
3.3.2
用正交函数族作最佳平方逼近
设 f ( x) C[a, b], span{0 ( x), 1 ( x), , n ( x)},
若 0 ( x), 1 ( x), , n ( x) 是满足条件(2.2)的正交函数族,
则
(i ( x), j ( x)) 0, i j
(k ( x), j ( x)) a j ( f ( x), k ( x)) (k 0,1,, n). (3.3)
n
而 j 0
b
j k.
(
),)
j (( xx))
(0x)dx 0,
j ((xx
( j , k )
)
j
k
a
j k.
Ak 0,
故法方程(3.3)的系数矩阵
(2.2)
Gn G(0 ( x), 1 ( x), , n ( x))
97
为非奇异对角阵,且方程(3.3)的解为
ak* ( f ( x), k ( x)) /( k (k ), k ( x))
(3.8)
(k 0,1,, n).
( k ( x), j ( x)) a j ( f ( x), k ( x))(k 0,1, , n). (3.3)
j 0
f ( x) C[a, b] 在 中的最佳平方逼近函数为
于是
n
n
S ( x)
*
k 0
( f ( x), k ( x))
k ( x)
2
k ( x).
(3.9)
2
98
由(3.5)可得均方误差为
n ( x)
2
f ( x) S n* ( x)
2
1
2
2
n
(
f
(
x
),
(
x
))
2
k
f ( x) 2
.
k ( x) 2
k 0
2
*
*
( x) 2 ( f ( x) S ( x), f ( x) S ( x))
(3.10)
由此可得贝塞尔(Bessel)不等式
( f ( x), f ( x)) ( S * ( x), f ( x))
n
(3.5)
*
2
n
(
a
(
x
)
)
f
(
x
)
. *
2
(3.11)
k
k
2
2
f ( x) 2 ak ( k ( x), f ( x)).
k 1
2
k 0
99
若 f ( x) C[a, b],按正交函数族 { k ( x)} 展开,系数
ak* (k 0,1, ) 按(3.8)计算,得级数
*
a
k k ( x)
(3.12)
k 0
系数 a k*
*
称这个级数为
的广义傅里叶(Foureir)级数,
f
(x
)
ak ( f ( x), k ( x)) /( k (k ), k ( x))
(3.8)
称为广义傅里叶系数. 它是傅里叶级数的直接推广.
(k 0,1,, n).
讨论特殊情况,设 {0 ( x), 1 ( x), , n ( x)} 是正交多
项式, span{0 ( x), 1 ( x), , n ( x)}, k ( x)( k 0,1,, n)
可由 1, x,, x n 正交化得到,则有下面的收敛定理.
100
定理8
设 f ( x) C[a, b], S * ( x) 是由(3.9)给出的
f (x)的最佳平方逼近多项式,其中 {k ( x), k 0,1,, n}
是正交多项式族,则有
lim f ( x) S n* ( x)
n
2
0.
( f ( x), k ( x))
S ( x)考虑函数
k ( x).按勒让德多项式
2
(3.9)
f ( (xx) )C[1, 1],
k 0
n
*
k
2
{P0 ( x), P1 ( x), , Pn ( x)} 展开,由(3.8),(3.9)可得
S n* ( x) a0* P0 ( x) a1* P1 ( x) an* Pn ( x),
(3.13)
101
n ( x)
其中
2
f ( x) S n* ( x)
2
1
2
P ( x))
( f ( x),
2nk (1f (1x), k ( x))
2
k
(3.1
a ( x)
..
f ( x) 2
f ( x) Pk ( x)dx
1 ( x )
( Pk ( x), Pk ( x)) k 02
k
2
(3.14)
2
*
k
根据均方误差公式(3.10),平方误差为
k ( x)
2
2
1
1
n
2
*2
f ( x)dx
ak .
k 0 2k 1
2
(3.15)
由定理8可得
lim f ( x) S n* ( x)
n
2
0.
102
*
*
* *一致收敛于
(3.13)
S n* ( xf) (x
)a满足光滑性条件,还有
P
(
x
)
a
P
(
x
)
a
如果
f (x)
0 0
1 1
S n n Pn ( x),
的结论.
定理9
设 f ( x) C 2 [1, 1], S n* ( x)由(3.13)给出,
则对任意 x [1, 1]和 0, 当 n充分大时有
f ( x) S n* ( x)
n
.
~ ,
公式(2.6)给出了首项系数为1的勒让德多项式 P
n
它具有以下性质.
~
Pn ( x)
n! d n
2
n
[(
x
1
)
].
n
(2n)! dx
103
(2.6)
定理10 在所有最高次项系数为1的 n次多项式中,
~
勒让德多项式 P
在 [1, 1] 上与零的平方误差最小.
n ( x)
证明 设 Qn (x)是任意一个最高次项系数为1的 n次
多项式,它可表示为
n 1
~
~
Qn ( x) Pn ( x) ak Pk ( x),
k 0
104
于是
Qn ( x)
2
2
1
(Qn ( x), Qn ( x)) Qn2 ( x)dx
1
n 1
~
~
~
2 ~
( Pn ( x), Pn ( x)) ak ( Pk ( x), Pk ( x))
k 0
~
~
( Pn ( x), Pn ( x))
~
Pn ( x)
2
2
当且仅当 a0 a1 an1 0 时等号才成立, 即当
~
Qn ( x) Pn ( x) 时平方误差最小.
105
例7
求 f ( x) e x 在 [1, 1]上的三次最佳平方逼近多
项式.
解
~
(
f
(
x
),
P
先计算
k ( x))
(k 0,1,2,3).
1
( f ( x), P0 ( x)) e dx e 2.3504;
1
e
1
x
1
( f ( x), P1 ( x)) xe x dx 2e 1 0.7358;
1
3
1
( f ( x), P2 ( x)) ( x 2 )e x dx
1 2
2
1
5 3 3
( f ( x), P3 ( x)) ( x x)e x dx
1 2
2
1
e
7
0.1431;
e
1
37 5e 0.02013.
e
106
由(3.14) 得
a0* ( f ( x), P0 ( x)) / 2 1.1752,
*
a
3x())f ( x),
P1( x
))1/ 2 1.1036,
( f ( x),1P
(
2
k
1
*
k
ak ( x )
f ( x) Pk ( x)dx.
1
( Pk ( x),*Pk ( x))
2
a2 5( f ( x), P2 ( x)) / 2 0.3578, (3.14)
a3* 7( f ( x), P3 ( x)) / 2 0.07046.
代入(3.13) 得三次最佳平方逼近多项式
S3* ( x) 0.9963 0.9979 x 0.5367 x 2 0.1761x 3 .
S n* ( x) a0* P0 ( x) a1* P1 ( x) an* Pn ( x),
(3.13)
107
均方误差
n ( x) 2 e S ( x)
x
*
3
2
3
2
*2
e dx
ak
1
k 0 2k 1
1
2x
0.0084.
最大误差
n ( x)
x
*
e
S
3 ( x)
0.0112.
如果 f ( x) C[a, b], 求 [ a, b]上的最佳平方逼近多项式,
做变换
ba
ba
x
t
2
2
(1 t 1),
108
于是
ba
ba
F (t ) f (
t
),
2
2
*
在 [1, 1]上可用勒让德多项式做最佳平方逼近多项式 S n (t ),
从而得到区间 [ a, b]上的最佳平方逼近多项式
S n* (
1
(2 x a b)).
ba
109
由于勒让德多项式 {Pk ( x)}是在区间 [1, 1]上由
{1, x,, x k ,} 正交化得到的,因此利用函数的勒让德展
开部分和得到最佳平方逼近多项式与由
S * ( x) a0 a1 x an x n
直接通过解法方程得到 H n 中的最佳平方逼近多项式是一致
的.
只是当 n 较大时法方程出现病态,计算误差较大,不能
使用,而用勒让德展开不用解线性方程组,不存在病态问题,
因此通常都用这种方法求最佳平方逼近多项式.
110
切比雪夫级数
3.3.3
如果 f ( x) C[1,1] ,按 {Tk ( x)}0 展成广义傅里叶级数,
由(3.12)可得级数
C0*
Ck*Tk ( x),
2
k 1
(3.16)
其中系数根据(3.8)式,由(2.12)得到
C
这里
*
k
2
1
f ( x)Tk ( x) dx
1
1 x2
,
Tk ( x) cos(k arccos x),
k 0,1, , (3.17)
x 1.
级数(3.16)称为 f (x) 在 [1,1] 上的切比雪夫级数.
111
若令 x cos , 0 ,则(3.16)就是 f (cos ) 的傅里
叶级数,其中
C
*
k
2
1
1
f (cos ) cos k d ,
k 0,1, .
(3.18)
根据傅里叶级数理论,只要 f (x) 在 [1,1]上分段连续,则
f (x)在 [1,1]上的切比雪夫级数(3.16)一致收敛于 f (x).
从而
C0*
f ( x)
Ck*Tk ( x).
2
k 1
(3.19)
112
取部分和
n
C0*
C ( x)
Ck*Tk ( x),
2
k 1
*
n
(3.20)
其误差为
f ( x) Cn* ( x) Cn*1Tn 1 ( x).
在 [1,1] 上 Tn1 ( x) 是均匀分布的,它的最大值 max Tn1 ( x)
1 x 1
最小,因此 Cn* ( x) 可作为 f (x) 在 [1,1] 上的近似最佳一致
逼近多项式.
113
例8
解
求 f ( x) e x 在 [1, 1]上的切比雪夫部分和 C3* ( x) .
由(3.18)得
C
*
k
2
0
e cos cos k d ,
k 0,1,2,3.
利用第4章的数值积分法得
C0* 2.53213176, C1* 1.13031821,
C2* 0.27149534, C3* 0.04433685.
由(3.20)及 Tk (x) 的表达式有
C3* ( x) 0.994531 0.997308 x 0.542991x 2 0.177347 x 3 ,
及 e x C3* ( x)
0.00607.
114
3.4
曲线拟合的最小二乘法
3.4.1
最小二乘法及其计算
在函数的最佳平方逼近中 f ( x) C[a, b], 如果 f (x)
只在一组离散点集 {xi , i 0,1, , m} 上给定,这就是科
学实验中经常见到的实验数据 {( xi , yi ), i 0,1,, m}的
曲线拟合.
115
问题为利用 yi f ( xi ), i 0,1,, m, 求出一个函数
y S * ( x) 与所给数据{( xi , yi ), i 0,1,, m} 拟合.
记误差
i S * ( xi ) yi ,
i 0,1,, m,
( 0 , 1 ,, m )T
116
设 0 ( x), 1 ( x), , n ( x)是 C[a, b] 上线性无关函数族,
在 span{0 ( x), 1 ( x), , n ( x)}中找一函数 S * ( x),
使误差平方和
2
2
m
m
[ S * ( xi ) yi ]2
i 0
min
2
i
S ( x )
i 0
m
2
[
S
(
x
)
y
]
i
i
(4.1)
i 0
这里
S ( x) a00 ( x) a11 ( x) ann ( x)
( n m)
(4.2)
117
这个问题称为最小二乘逼近,几何上称为曲线拟合的
最小二乘法.
用最小二乘求拟合曲线时,首先要确定 S (x ) 的形式.
确定 S (x ) 的形式问题不仅是数学问题, 还与问题的
实际背景有关.
通常要用问题的运动规律及给定的数据进行数据描图,
确定 S (x )的形式, 然后通过实际计算选出较好的结果.
118
S (x ) 的一般表达式为(4.2)表示的线性形式.
若 k (x)是 k 次多项式,S (x ) 就是 n 次多项式.
为了使问题的提法更有一般性,通常在最小二乘法中
S ( x) a00 ( x) a11 ( x) ann ( x) (n m)
考虑加权平方和
(4.2)
2
2
m
( xi )[ S ( xi ) yi ]2 .
(4.3)
i 0
这里 ( x) 0是 [ a, b]上的权函数,它表示不同点 ( xi , f ( xi ))
处的数据比重不同.
119
用最小二乘法求拟合曲线的问题,就是在形如(4.2)的
S (x )中求一函数 y S * ( x),使(4.3)取得最小.
这样,最小二乘问题就转化为求多元函数
( n m)
m
n
(4.2)
( xi )[ a j j ( xi ) f ( xi )]2
(4.4)
S (Ix()a
a1
) n)a11 ( x) ann ( x)
,
,a
0 (x
0, 0
i 0
j 0
m
*
2
,
,
a
的极小点 (a , a
n()x问题.
i )[ S ( xi ) yi ] .
*
0
2*
1
2
(4.3)
i 0
由求多元函数极值的必要条件,有
m
n
I
2 ( xi )[ a j j ( xi ) f ( xi )] k ( xi ) 0
ak
i 0
j 0
(k 0,1,, n).
120
若记
m
( j , k ) ( xi ) j ( xi ) k ( xi ),
(4.5)
i 0
m
( f , k ) ( xi ) f ( xi ) k ( xi ) d k
i 0
(k 0,1, , n).
上式可改写为
n
(
j 0
k
, j )a j d k
(k 0,1,, n). (4.6)
这方程称为法方程,可写成矩阵形式
121
Ga d .
T
T
a
(
a
,
a
,
,
a
)
,
d
(
d
,
d
,
,
d
)
,
其中
0
1
n
0
1
n
( 0 , 0 )
( , )
G 1 0
( n , 0 )
( 0 , 1 )
(1 , 1 )
( n , 1 )
( 0 , n )
(1 , n )
.
( n , n )
(4.7)
要使法方程(4.6)有唯一解,就要求矩阵 G非奇异,
而 0 ( x), 1 ( x), , n ( x)在 [ a, b]上线性无关不能推出
n
矩阵 G非奇异,必须加上另外的条件.
(k , j )a j d k (k 0,1,, n). (4.6)
j 0
122
定义7
设 0 ( x), 1 ( x), , n ( x) [a, b]的任意线
性组合在点集 {xi , i 0,1,, m}(m n) 上至多只有 n 个
不同的零点,则称 0 ( x), 1 ( x), , n ( x)在点集
{xi , i 0,1,, m}上满足哈尔(Haar)条件.
显然 1, x,, x n在任意 m(m n)个点上满足哈尔条件.
如果 0 ( x), 1 ( x), , n ( x) [a, b] 在 {xi }0m 上满足
哈尔条件,则法方程(4.6) 的系数矩阵(4.7) 非奇异,于是
方程(5.6)存在唯一的解 ak ak* , k 0,1, , n. 从而得到
n
函数 f (x) 的最小二乘解为
(k , j )a j d k
j 0
(k 0,1,, n). (4.6)
123
S * ( x) a0*0 ( x) a1*1 ( x) an* n ( x).
这样得到的 S * ( x), 对任何形如(4.2)的 S (x ) ,都有
m
m
i 0
i 0
*
2
2
(
x
)[
S
(
x
)
f
(
x
)]
(
x
)[
S
(
x
)
f
(
x
)]
,
i
i
i
i
i
i
故 S * ( x)确是所求最小二乘解.
S ( x) a00 ( x) a11 ( x) ann ( x)
( n m)
(4.2)
124
给定 f (x) 的离散数据 {( xi , yi ), i 0,1,, m},
一般可取 span{1, x,, x n } ,但这样做当 n 3时,
求解法方程(4.6)将出现系数矩阵 G为病态的问题,
通常对 n 1的简单情形都可通过求法方程(4.6)得到
S * ( x).
有时根据给定数据图形,其拟合函数 y f (x) 表面上
不是(4.2)的形式,但通过变换仍可化为线性模型.
例如, S ( x) aebx ,若两边取对数得
S ( x) a00 ( x) a11 ( x) ann ( x)
ln S ( x) ln a bx,
( n m)
125
(4.2)
此时,若令 S ( x) ln S ( x),
A ln a,
B b,
则
S ( x) A Bx ,
这样就变成了形如(4.2)的线性模型 .
126
例9
解
已知一组实验数据如下,求它的拟合曲线.
xi
1
2
3
4
5
fi
4
4.5
6
8
8.5
i
2
1
3
1
1
将所给数据在坐标纸上标出,见图3-5.
127
图3-5
从图中看到各点在一条直线附近,故可选择线性函数
作拟合曲线,
令 S1 ( x) a0 a1 x, 这里 m 4, n 1, 0 ( x) 1,
1 ( x) x, 故
4
( 0 , 0 ) i 8,
i 0
4
( 0 , 1 ) (1 , 0 ) i xi 22,
i 0
4
(1 , 1 ) i xi2 74,
i 0
4
( 0 , f ) i f i 47,
i 0
4
(1 , f ) i xi f i 145.5.
i 0
128
由(4.6)得方程组
8a0 22a1 47
22a0 74a1 145.5.
n
(
, j )a j d k
(k 0,1,, n). (4.6)
j 0
解得 a0 2.77, a1 1.13. 于是所求拟合曲线为
k
S1* ( x) 2.77 1.13x.
129
关于多项式拟合,Matlab中有现成的程序
a polyfit ( x, y, m)
其中输入参数 x, y为要拟合的数据, m为拟合多项式的次数,
输出参数 a 为拟合多项式的系数.
利用下面的程序,可在Matlab中完成上例的多项式拟合.
130
x=[1 2 3 4 5];
f=[4 4.5 6 8 8.5];
aa=polyfit(x,f,1);
y=polyval(aa,x);
plot(x,f,’r+’,x,y,’k’)
xlabel(‘x’);
ylabel(‘y’);
gtext(‘y=s1(x)’)
131
结果如下:
132
例10 设数据 ( xi , yi )(i 0,1,2,3,4) 由表3-2给出,
表中第4行为 ln yi yi ,通过描点可以看出数学模型为
y aebx , 用最小二乘法确定 a及 b .
表3 2
i
0
1
2
3
4
xi
1.00
1.25
1.50
1.75
2.00
yi
5.10
5.79
6.53
7.45
8.46
133
解
用给定数据描图可确定拟合曲线方程为 y aebx ,
它不是线性形式. 两边取对数得
ln y ln a bx,
若令 y ln y, A ln a, 则得 y A bx, {1, x}.
为确定 A, b ,先将 ( xi , yi )转化为 ( xi , yi ), 数据表见表3-2.
0 ( x) 1, 1 ( x) x, ( x) 1,
得
表3 2
(0 , 0 ) 5,
41
i
0
2
( 0 , 1 ) xi 7.5,
xi
1.00
1.25
1.50
i 0
3
4
1.75
2.00
5.10
54 .792
6.53
7.45
(1 , 1 ) xi 11.875,
yi 1.629 1i .0756 1.876 2.008
yi
8.46
2.135
134
4
( 0 , y ) yi 9.404,
i 0
4
(1 , y ) xi yi 14.422.
i 0
故有法方程
5 A 7.50b 9.404
7.50 A 11.875b 14.422.
解得
A 1.122, b 0.505, a e A 3.071.
于是得最小二乘拟合曲线为
y 3.071e0.505x .
135
利用下面的程序,可在Matlab中完成曲线拟合.
x=[1.00 1.25 1.50 1.75 2.00];
y=[5.10 5.79 6.53 7.45 8.46];
y1=log(y);
aa=polyfit(x,y1,1);
a=aa(1); b=exp(aa(2));
y2=b*exp(a*x);
plot(x,y,’r+’,x,y2,’k’)
xlabel(‘x’);
ylabel(‘y’);
gtext(‘y=a*exp(bx))’;
136
结果如下:
137
3.4.2
用正交多项式做最小二乘拟合
用最小二乘法得到的法方程组(4.6),其系数矩阵
G是病态的.
如果
0 ( x), 1 ( x), , n ( x) 是关于点集
n
( k , j )a j d k
(k 0,1,, n). (4.6)
{xi }j 0
(i 0,1, , m) 带权 ( xi ) (i 0,1,, m)
函数族,即
0,
( j , k ) ( xi ) j ( xi ) k ( xi )
i 0
Ak 0,
m
j k,
j k,
(4.8)
138
则方程(4.6)的解为
m
ak* n
j 0
( x ) f ( x )
( xi )
( f ,k )
i 0 m
(k 0,1,, n).
k,,
k ))a d
((
k
j
j
k ( x(i k)k2 (0x,i1),, n). (4.6)
i
i
k
i 0
(4.9)
且平方误差为
2
2
f
2
2
n
Ak (ak* ) 2 .
k 0
139
接下来根据给定节点 x0 , x1 , , xm 及权函数 ( x) 0,
构造带权 ( x ) 正交的多项式 {Pn ( x)} .
注意 n m,用递推公式表示 Pk (x) ,即
P0 ( x) 1,
(4.10)
P1 ( x) ( x 1 ) P0 ( x),
P ( x) ( x ) P ( x) P ( x)
k 1
k
k k 1
k 1
(k 1,2,, n 1).
这里 Pk (x) 是首项系数为1的 k 次多项式, 根据 Pk (x) 的
正交性,得
140
m
2
(
x
)
x
P
i i k ( xi )
( xPk ( x), Pk ( x))
k 1 i 0m
( Pk ( x ), Pk ( x ))
2
(
x
)
P
(
x
)
i
k
i
i 0
( xPk , Pk )
( Pk , Pk )
m
( xi ) Pk2 ( xi )
( Pk , Pk )
i 0
k
m
( Pk 1 , Pk 1 )
2
(
x
)
P
(
x
)
i k 1 i
i 0
(4.11)
(k 1,2,, n 1).
下面用归纳法证明这样给出的 {Pk ( x)}是正交的.
141
由(4.10)第二式及(4.11)中 1 的表达式,有
( P0 , P1 ) ( P0 , xP0 ) 1 ( P0 , P0 )
( P0 , xP0 )
( xP0 , P0 )
( P0 , P0(4.10)
) 0.
( P0 , P0 )
P0 ( x) 1,
P1 ( x) ( x 1 ) P0 ( x),
( Pl , Ps ) 0(l s) 对 s 0,1,, l 1 及 l 0,1, , k ;
假定
Pk 1 ( x) ( x k 1 ) Pk ( x) k Pk 1 ( x)
k n 均成立,要证 ( Pk 1 , Ps ) 0 对 s 0,1, , k 均成立.
(k 1,2,, n 1).
由(4.10)有
( Pk 1 , Ps ) (( x k 1 ) Pk , Ps ) k ( Pk 1 , Ps )
(4.12)
P0 ( x) 1,
(4.10)
0
x k ,P1s))P
( xk),
1 ( Pk , Ps ) k ( Pk 1 , Ps ).
P1 ( x) ( xP
P ( x) ( x ) P ( x) P ( x)
142
k 1
k
k k 1
k 1
(k 1,2,, n 1).
由归纳法假定,当 0 s k 2 时
( Pl , Ps ) 0,
( Pk 1 , Ps ) 0.
另外,xPs (x) 是首项系数为1的 s 1 次多项式,它可由
P0 , P1 , , Ps 1 的线性组合表示.
而 s 1 k 1,由归纳法假定又有
( xPk , Ps ) ( Pk , xPs ) 0,
于是由(4.12),当 s k 2 时,( Pk 1 , Ps ) 0.
( Pk 1 , Ps ) (( x k 1 ) Pk , Ps ) k ( Pk 1 , Ps )
(4.12)
( xPk , Ps ) k 1 ( Pk , Ps ) k ( Pk 1 , Ps ).
143
再考虑
( Pk 1 , Pk 1 ) ( xPk , Pk 1 ) k 1 ( Pk , Pk 1 ) k ( Pk 1 , Pk 1 ),
(4.13)
由假定有
( Pk , Pk 1 ) 0,
k 1
( xPk , Pk 1 ) ( Pk , xPk 1 ) ( Pk , Pk c j Pj ) ( Pk , Pk ).
j 0
利用(4.11)中 k 表达式及以上结果,得
( Pk 1 , Pk 1 ) ( xPk , Pk 1 ) k ( Pk 1 , Pk 1 )
( Pk , Pk )
k
).
( Pk , Pk ) ( P(kP
,kP
1
k 10
k ,)P
144
最后,由 k 和 k 的表达式(4.11)有
( Pk 1 , Pk ) ( xPk , Pk ) k 1 ( Pk , Pk ) k ( Pk , Pk 1 )
( xPk ,(P
xP
P)
k)
k,
( xPk , Pk )
( Pkk , Pk ) 0.
k 1
( Pk , P
(kP)k , Pk )
至此已证明了由(4.10)及(4.11)确定的多项式
( Pk , Pk )
{Pk ( x)}
k
( Pk 1 , Pk 1 )
{
x
}
组成一个关于点集 i 的正交系.
用正交多项式 {Pk ( x)}的线性组合作最小二乘曲线拟合,
只要根据公式(4.10)及(4.11)逐步求 Pk (x) 的同时,
*
相应计算出系数 a k
145
m
( f , Pk )
a
( Pk , Pk )
*
k
(x ) f (x )P (x )
i
i 0
i
k
m
( x )P
i 0
i
2
k
i
(k 0,1, , n),
( xi )
并逐步把 ak* Pk ( x)累加到 S (x ) 中去,最后就可得到所求的
拟合曲线
y S ( x) a0* P0 ( x) a1* P1 ( x) an* Pn ( x).
这里 n可事先给定或在计算过程中根据误差确定.
用这种方法编程序不用解方程组,只用递推公式,并
且当逼近次数增加一次时,只要把程序中循环数加1,其余
不用改变.
146
147
2020/5/1
147
2020/5/1
148
148
149
2020/5/1
149
150
151
3.5
有 理 逼 近
3.5.1
有理逼近与连分式
有理函数逼近是指用形如
n
Rnm ( x)
的函数逼近 f (x ).
Pn ( x)
Qm ( x)
a
k
xk
k
xk
b
k 0
与前面讨论一样,如果 f ( x) Rnm ( x)
最佳有理一致逼近.
(5.1)
k 0
m
最小就可得到
152
如果 f ( x) Rnm ( x)
2
最小则可得到最佳有理平方逼近
函数.
本节主要讨论利用函数的泰勒展开获得有理逼近函数
的方法.
对函数 ln( 1 x)用泰勒展开得
ln( 1 x) (1)
k 1
k 1
xk
k
x [1,1].
(5.2)
取部分和
n
S n ( x) (1) k 1
k 1
xk
ln( 1 x).
k
153
另一方面若对(5.2)式用辗转相除可得到 ln( 1 x) 的
一种连分式展开
x
1 x
1
k
1 x
x
k
1
2
ln( 1 x) (1)
x2
2[x1,1].
3k
k 1
22 x
4
5
ln( 1 x )
(5.2)
x 1 x 1 x 2 2 x 22 x
.
1 2 3 4 5
(5.3)
154
(5.3)右端为 ln( 1 x) 的无穷连分式的前5项,最后式子
是它的紧凑形式.
若取(5.3)的前2,4,6,8项,则可分别得到 ln( 1 x)
的以下有理逼近
2x
6 x 3x 2
R22 ( x)
,
R11 ( x ) 2 x ,
2
6 6x x
60 x 60 x 2 11x 3
R33 ( x)
,
2
3
60 90 x 36 x 3 x
2 260 x 3 25 x 4
420
x
630
x
R44 ( x)
. (5.4)
2
3
4
420 840 x 540 x 120 x 6 x
155
若用同样多项的泰勒展开部分和 S 2 n ( x) 逼近 ln( 1 x)
并计算 x 1 处的值 S 2n (1)及 Rnn (1),计算结果见表3-3.
表3 3
n
S 2 n (1)
s ln( 2) S 2 n (1)
Rnn (1)
R ln( 2) Rnn (1)
1
0.50
0.19
0.667
0.026
2
0.58
0.11
0.69231
0.00084
3
0.617
0.076
0.693122
0.000025
4
0.634
0.058
0.69314642
0.00000076
ln 2 的准确值为 0.69314718, 从表3-1可以看出,
R44 (1) 0.69314642,
S8 (1) 0.634,
156
由此看出 R44 (1) 的精度比 S8 (1) 高出近10万倍,
但它们的计算量是相当的,这说明用有理逼近比多项
式逼近好得多.
157
例11
给出有理函数
2 x 4 45 x 3 381x 2 1353x 1511
R43 ( x)
x 3 21x 2 157 x 409
用辗转相除法将它化为连分式并写成紧凑形式.
解
用辗转相除可逐步得到
158
4 x 2 64 x 284
R43 ( x) 2 x 3
x 3 21x 2 157 x 409
2x 3
4
6( x 9)
x5
x 2 16 x 71
4
2x 3
x5
6
x7
2x 3
8
x9
4
6
8
.
x5 x7 x9
本例中用连分式计算 R43 ( x)的值只需3次除法,1次乘
法和7次加法.
159
若直接用多项式计算的秦九韶算法则需6次乘法和1次
除法及7次加法.
可见将 Rnm (x) 化成连分式可节省计算乘除法次数.
对一般的有理函数,(5.1)可转化为一个连分式
cl
c2
Rnm ( x) P1 ( x)
.
x d1 x d l
它的乘除法运算只需 max( m, n) 次.
而直接用有理函数(5.1)计算乘除法次数为 n m 次.
160
3.5.2
帕德逼近
利用函数 f (x) 的泰勒展开可以得到它的有理逼近.
设 f (x)在 x 0的泰勒展开为
N
1
f ( x)
f
k 0 k!
(k )
f ( N 1) ( ) N 1
(0) x
x
.
( N 1)!
k
(5.5)
它的部分和记作
N
1
P( x)
f
k 0 k!
(k )
(0) x k
N
c
k
xk .
(5.6)
k 0
161
定义8
设 f ( x) C N 1 (a, a), N n m,
如果有理函数
a0 a1 x an x n
Pn ( x)
Rnm ( x)
,
m
1 b1 x bm x
Qm ( x)
(5.7)
其中 Pn ( x), Qm ( x) 无公因式,且满足条件
(k )
Rnm
(0) f
(k )
(0)
(k 0,1, , N ),
(5.8)
则称 Rnm (x) 为函数 f (x) 在 x 0 处的 (n, m) 阶帕德逼近,
记作 R(n, m) ,简称 R(n, m) 的帕德逼近.
162
根据定义,若令
h( x) P( x)Qm ( x) Pn ( x),
则满足条件(5.8)等价于
h ( k ) (0) 0
即
k 0,1,, N .
(k )
(5.8)
Rnm
(0) f ( k ) (0)
(k 0,1, , N ),
(k )
(k )
h (0) ( P( x)Qm ( x) Pn ( x))
0,
k 0,1,, N .
x 0
由于 Pn( k ) (0) k!ak , 应用莱布尼兹求导公式得
( P( x)Qm ( x) Pn ( x)) ( k )
k
x 0
k! c j bk j k!ak 0
j 0
k 0,1, , N , 163
1
这里 c j f
j!
( j)
(0) 是由(5.6)得到的,上式两端除 k! ,
并由 b0 1, b j 0(当j m时), 可得
k 1
ak c j bk j ck
(k 0,1, , n)
(5.9)
j 0
N
及
1
P( x)
f
k 0 k!
k 1
c j bk j ck
(k )
(0) x k
N
c
k
xk .
(5.6)
k 0
(k n 1, , n m)
(5.10)
j 0
注意当 j m 时 b j 0, 故(5.10)可写成
164
cn m 1bm cn 1b2 cn b1 cn 1 ,
c
n m 2 bm cn b2 cn 1b1 cn 2 ,
cn bm cn m 2b2 cn m 1b1 cn m ,
(5.11)
, 若记
其中 j 0 时 c j 0,
cn m 1
c
H nm 2
cn
cn 1
cn
cn m 2
b (bm , bm1 ,, b1 )T ,
cn
cn 1
, (5.12)
cn m 1
c (cn1 , cn 2 ,, cn m )T .
165
则方程组(5.11)的矩阵形式为
Hb c .
定理11
设 f ( x) C N 1 (a, a),
N n m, 则形如
(5.7)的有理函数 Rnm (x) 是 f (x) 的 (n, m)阶帕德逼近的
充分必要条件是多项式 Pn ( x), Qm ( x) 的系数 a0 , a1 ,, an
及 b0 , b1 ,, bm 满足方程组(5.9)及(5.11).
166
根据定理11, 求 f (x) 的帕德逼近时,
首先要由(5.11)解出 Qm (x)的系数 b0 , b1 ,, bm ,
再由(5.9)直接算出 Pn (x) 的系数 a0 , a1 ,, an .
f (x)的各阶帕德逼近可列成一张表,称为帕德表(见表
cn m 1bm cn 1b2 cn b1 cn 1 ,
c
3-4).
n m 2 bm cn b2 cn 1b1 cn 2 ,
表3 - 4
m
0
1
2
3
n cn bm cn m 2b2 cn m 1b1 cn m ,
(5.11)
4
0
(0,0)
(0,1)
(0, 2)
(0,3)
(0, 4)
1
(1,0)
(1,1)
(1, 2)
(1,3)
(1, 4)
( 2,0)
( 2,1)
( 2, 2)
( 2,3)
( 2, 4)
( 4,0)
( 4,1)
( 4, 2)
( 4,3)
( 4, 4)
2
k 1
ak 3 c(3
,2
0,)1, (,3n,)3)
jb
)
j c(k3,1) ( k(3
,k0
4
j 0
((5.9)
3, 4)
167
例12 求 f ( x) ln( 1 x)的帕德逼近 R (2,2) 及 R(3,3).
解
由 ln( 1 x) 的泰勒展开
1 2 1 3 1 4
ln( 1 x) x x x x
2
3
4
1
1
1
得 c0 0, c1 1, c2 , c3 , c4 , .
2
3
4
当 n m 2时,由(5.11)得
1
1
b
b
,
2
1
2
3
1
1
1
b2 b1 .
3
4
2
求得 b1 1, b2 1 , 再由(5.9)得
6
168
a0 0,
a1 1,
a2
1
,
2
于是得
1 2
x
6 x 3x 2
2
.
R22 ( x)
2
1
6 6x x
1 x x2
6
x
当 n m 3时,由(5.11)得
1
1
1
b
b
b
,
3
2
1
2
3
4
1
1
1
1
b
b
b
,
3
2
1
3
4
5
2
1 b 1 b 1 b 1 ,
3
2
1
4
5
6
3
169
解得
b1
3
,
2
b2
3
,
5
b3
1
.
20
代入(5.9)得
a0 0,
a1 1,
a2 1,
11
a3
.
60
于是得
11 3
xx
x
60 x 60 x 2 11x 3
60
R33 ( x )
.
2
3
3
3 2
1 3
60 90 x 36 x 3 x
1 x x
x
2
5
20
2
170
可以看到这里得到的 R22 ( x) 及 R33 ( x) 与 ln( 1 x) 的前面
连分式展开得到的有理逼近(5.4)结果一样.
为了求帕德逼近 Rnm (x) 的误差估计,由(5.9)及(5.11)
求得的 Pn ( x), Qm ( x)系数 a0 , a1 ,, an 及 b0 , b1 ,, bm ,直
接代入则得
f ( x)Qm ( x) Pn ( x) x
n m 1
m
( bk cn m 1l k )x l ,
l 0 k 0
将 Qm (x)除上式两端,即得
171
x
n m 1
f ( x) Rnm ( x)
l
r
x
l
l 0
(5.13)
,
Qm ( x)
m
其中 rl bk cn m l 1 k .
k 0
当 x 1时可得误差近似表达式
f ( x) Rnm ( x) r0 x
n m 1
m
,
r0 bk cn m 1 k .
k 0
172
3.6
三角多项式逼近与快速傅里叶变换
当模型数据具有周期性时,用三角函数特别是正弦和余
弦函数作为基函数是合适的,这时前面讨论的用多项式、分
段多项式或有理函数作基函数都是不适合的.
用正弦函数和余弦函数级数表示函数称为傅里叶变换
(简称傅氏变换).
173
3.6.1
最佳平方三角逼近与三角插值
设 f (x) 是以 2π
1
S n ( x) a0 a1 cos x b1 sin x an cos nx bn sin nx
2
(6.1)
做最佳平方逼近函数. 由于三角函数族
1,cos x,sin x, ,cos kx,sin kx,
在 [0,2 π] 上是正交函数族,于是 f (x) 在 [0,2 π] 上的最小
平方三角逼近多项式 S n (x)
174
1
ak
π
1
bk
π
2
0
f ( x) cos kxdx
(k 0,1, , n),
(6.2)
2
0
f ( x) sin kxdx
(k 0,1, , n),
ak , bk 称为傅里叶系数.
函数 f (x) 按傅里叶系数展开得到的级数
1
a0 (ak cos kx bk sin kx)
2
k 1
(6.3)
就称为傅里叶级数.
175
只要 f (x )在 [0,2 π] 上分段连续,则级数(6.3)一致
收敛到 f (x).
对于最佳平方逼近多项式(6.1)有
1
a0 (ak cos2 kx bk sin
kx)
2
2 f ( x)k
1 S ( x )
f ( x) S
n
2
2
2
n
(6.3)
( x) 2 .
1
由此可以得到相应于(3.11)的贝塞尔不等式
S n ( x) a0 a1 cos x b1 sin x an cos nx bn sin nx
2
n
(6.1)
1 2
1 2
2
2
2
a0 (ak bk )
[ f ( x)] dx.
2
*
2
(
a
(
x
)
)
f ( x) 2 .
k
k
2
因为右边不依赖于
k 1
n ,左边单调有界,所以级数
n
2
k 1
0
(3.11)
176
1 2
a0 ( ak2 bk2 )
2
k 1
收敛,并有
lim ak lim bk 0.
k
k
当 f (x) 只在给定的离散点集
2π
x
j
,
j
0
,
1
,
,
N
1
j
N
上已知时,则可类似得到离散点集正交性与相应的离散傅
里叶系数.
下面只给出奇数个点的情形.
177
令
xj
2 πj
2m 1
( j 0,1,,2m),
可以证明对任何 0 k , l m 成立
2m
0,
sin l x j sin kx j 2m 1
,
j 0
2
0,
2m 1
2m
cos l x j cos kx j
j 0
2
2m 1,
2m
0 k,
cos l x j sin kx j 0,
j 0
l k , l k 0;
l k 0;
l k;
l k 0;
l k 0;
j m.
178
这表明函数族 {1, cos x,sin x, ,cos mx,sin mx} 在点集
{x j
2 πj
}
2m 1
上正交.
若令 f j f ( x j )
( j 0,1,,2m), 则 f (x) 的最小二
乘三角逼近为
n
1
S n ( x) a0 (ak cos kx bk sin kx),
2
k 1
n m,
其中
2m
2
2πjk
ak
f
cos
j
2m 1 j 0
2m 1
(k 0,1, , m),
179
2m
2
2 πjk
bk
f
sin
j 2m 1
2m 1 j 0
(k 1, , n).
(6.4)
当 n m时
Sm ( x j ) f j
( j 0,1, ,2m) ,
于是
m
1
S m ( x) a0 (ak cos kx bk sin kx)
2
k 1
就是三角插值多项式,系数仍由(6.4)表示.
180
一般情形,假定 f (x) 是以 2π为周期的复函数,给定
在 N个等分点 x j
2π
2π
j ( j 0,1, , N 1) 上的值 f j f
N
N
j ,
由于
eijx cos( jx) i sin( jx) ( j 0,1,, N 1, i 1),
所以函数族 {1, eix ,, ei ( N 1) x }在区间 [0,2 π] 上是正交的.
函数 e ijx 在等距点集 xk 2π k (k 0,1, , N 1) 上的值
N
e ijxk 组成的向量记作
j (1, e
ij
2π
N
,, e
ij
2π
( N 1)
N
)T .
181
当 j 0,1,, N 1 时,N 个复向量 0 , 1 , , N 1 具有
如下正交性:
N 1
(l , s ) e
il
2π
k
N
e
is
2π
k
N
k 0
N 1
e
i(l s )
2π
k
N
k 0
0, l s;
N , l s.
(6.5)
182
事实上,令 r e
i (l s )
0 l N 1,
2π
N
, 若 0 l , s N 1, 则有
( N 1) s 0,
于是
( N 1) l s N 1,
即
N 1 l s
N 1
1
1;
N
N
N
若 l s 0, 则 r 1, 从而
r N ei ( l s ) 2 π 1;
183
于是
(l , s )
N 1
r
k 0
k
1 r N
1 r
0.
若 l s, 则 r 1, 于是
(s , s )
N 1
k
r
N.
k 0
这就证明了(6.5)成立. 即 0 , 1 , , N 1 是正交的.
2π
因此, f (x)在 N个点 {x j
j ( j 0,1, , N 1)} 上的
N
0, l s;
最小二乘傅里叶逼近为
(l ,s )
(6.5)
N
,
l
s
.
184
S ( x)
n 1
ikx
c
e
k ,
n N,
(6.6)
k 0
其中
1
ck
N
N 1
f
j 0
j
e
-ikj
2π
N
,
(k 0,1, , n 1). (6.7)
在(6.6)中,若 n N ,则 S (x )为 f (x) 在点
x j ( j 0,1, , N 1) 上的插值函数, 即 S ( x j ) f ( x j ),
于是由(6.6)得
N 1
f j ck e
k 0
ik
2π
j
N
,
( j 0,1, , N 1).
(6.8)
185
(6.7)是由 { f j } 求 {ck } 的过程,称为 f (x) 的离散
傅里叶变换. 简称DFT,
而(6.8)是由 {ck }求 { f j } 的过程,称为反变换.
1
ck
N
N 1
f
j 0
N 1
f j ck e
j
ik
e
-ikj
2π
j
N
2π
N
,
,
(k 0,1, , n 1). (6.7)
( j 0,1, , N 1).
(6.8)
k 0
186
3.6.2
N 点DFT与FFT算法
不论是按(6.7)式由 { f j } 求 {ck } ,或是按(6.8)
由 {ck } 求 { f j },还是由(6.4)计算傅里叶逼近系数 ak , bk ,
都可归结为计算
N 1
2 kj 2 m
2 πjk
c j ak xk N
, f (j cos
j 0,1, , N 1),
2m 1 (6.4)
k 0 2m 1 j 0
(6.9)
2m
其中 {xk }0N 1为已知的输入数据,
2
2{
πcjkj }0N 1为输出数据,而
bk
N e
i 2N
f
2m 1
j
sin
2m 1
2j0
2
cos
isin
, i 1.
N
N
187
(6.9)式称为 N点DFT,表面上看计算 c j 需要 N次复
数乘法和 N次复数加法,称为 N个操作,计算全部 c j 共要 N 2
个操作.
N 1
(6.9)
c 当 Nx 较大且处理数据很多时,就是用高速的电子计算
wkj , ( j 0,1, , N 1),
j
k 0
k
机,很多实际问题仍然无法计算,
直到20世纪60年代中期产生了FFT算法,大大提高了
运算速度,才使DFT得到更广泛的应用.
FFT是快速算法的一个典范,其基本思想就是尽量减
少乘法次数.
188
对于任意正整数 k, j ,成立
Nj Nk Nj k , NjN k Nk (周期性),
jk
Njk N / 2 Njk (对称性), jN
Nk .
由周期性可知所有 Njk ( j, k 0,1,, N 1) 中,做多有 N个不
同的值 N0 , 1N ,, NN 1.
特别地,有
N0 NN 1, NN / 2 1.
当 N 2 p 时, Njk只有 N / 2个不同的值.
利用这些性质可将(6.9)式对半折成两个和式,再将
对应项相加,有
cj
N / 2 1
x
k 0
k
jk
N
N / 2 1
x
k 0
N / 2 k
j ( N / 2 k )
N
N / 2 1
[ x
k 0
k
(1) j xN / 2 k ] Njk .189
依下标奇、偶分别考察,则
c2 j
N / 2 1
(x
c2 j 1
若令
k 0
k
x N / 2 k ) Njk/ 2 ,
N / 2 1
jk
(
x
x
)
k N / 2 k N / 2 .
k 0
yk ( xk xN / 2 k ),
y N / 2 k ( xk xN / 2 k ) Nk ,
则可将 N 点DFT归结为两个 N / 2点的DFT:
N / 2 1
jk
c
y
2
j
k
N
/2,
k 0
N / 2 1
jk
c
y
2 j 1
N / 2 k N / 2 .
k 0
j 0,1, , N / 2 1,
如此反复施行二分手续即可得到FFT算法.
190
以
N 23 为例说明FFT的计算方法,此时 k , j 0,1,, N 1 7
在(6.9)中将 N 8 记为 ,于是(6.9)式的和为
cj
7
x
k 0
k
jk
,
( j 0,1, ,7).
(6.10)
将 k, j 用二进制表示为
k k 2 2 2 k1 21 k0 20 (k 2 k1k0 ),
j j2 2 2 j1 21 j0 20 ( j2 j1 j0 );
其中 kr , jr (r 0,1,2) 只能取0或1.
例如 6 22 22 0 20 (110).
根据 k, j 表示法,有
c j c( j2 j1 j0 ),
xk x(k 2 k1k0 ).
191
公式(6.10)可表示为
7
c j xk kj , (6.10)
2
1
0
c( j2 j1 j0 ) x(k 2kk10k0 ) ( k 2 k1k0 )( j2 2 j1 2 j0 2 )
1
1
1
(6.11)
k 0 0 k1 0 k 2 0
j1 ( k1k0 0 )
1 1
j0 ( k 2 k1k 0 )
j ( k 00)
x(k 2 k1k0 )
2 0
k0 0
k1 0 k 2 0
1
若引入记号
A0 (k2 k1k0 ) x(k2 k1k0 ),
1
A1 (k1k0 j0 ) A0 (k 2 k1k0 ) w j0 ( k 2 k1k0 ) ,
k2 0
1
j1 ( k1k 0 0 )
A
(
k
j
j
)
A
(
k
k
j
)
w
,
1
1 0 0
2 0 1 0
k1 0
1
j 2 ( k 0 00)
A ( j j j )
A
(
k
j
j
)
w
.
3
2 1 0
2
0 1 0
k0 0
(6.12)
192
则(6.11)变成
c( j2 j1 j0 ) A3 ( j2 j1 j0 ).
若注意 j 2 j N / 2 (1) j , 公式(6.12)还可进一步
p1
0
0
0
简化为
A1 (k1k0 j0 )
1
A (k
k 2 0
0
2
k1k0 ) j0 ( k 2 k1k0 )
A0 (0k1k0 )
j0 ( 0 k1k 0 )
A0 (1k1k0 )
j0 2 2
j
0 ( 0 k1k 0 )
[ A0 (0k1k0 ) (1) j0 A0 (1k1k0 )] j0 ( 0 k1k0 ) ,
A1 (k1k0 0) A0 (0k1k0 ) A0 (1k1k0 ),
A1 (k1k01) [ A0 (0k1k0 ) A0 (1k1k0 )] ( 0 k1k0 ) .
193
将这表达式中二进制表示还原为十进制表示:
k (0k1k0 ) k1 21 k0 20 ,
即 k 0,1,2,3, 得
2
A1 (2k ) A0 (k ) A0 (k 2 ),
2
k
A1 ( 2k 1) [ A0 ( k ) A0 ( k 2 )]w
( k 0,1,2,3).
(6.13)
同样(6.12)中的 A2 也可简化为
A2 (k0 j1 j0 ) [ A1 (0k0 j0 ) (1) j1 A1 (1k0 j0 )] j1 ( 0 k0 0 ) ,
即
A2 (k0 0 j0 ) A1 (0k0 j0 ) A1 (1k0 j0 ),
A2 (k01 j0 ) [ A1 (0k0 j0 ) A0 (1k0 j0 )] ( 0 k0 0 ) .
194
把二进制表示还原为十进制表示,得
A2 (k 22 j ) A1 (2k j ) A1 (2k j 22 ),
2
2
2k
A2 ( k 2 j 2) [ A1 ( 2k j ) A1 ( 2k j 2 )]w
(6.14)
( k 0,1; j 0,1).
同理(6.12)中 A3 可简化为
A3 ( j2 j1 j0 ) A2 (0 j1 j0 ) (1) j2 A2 (1 j1 j0 ),
即
A3 (0 j1 j0 ) A2 (0 j1 j0 ) A2 (1 j1 j0 ),
A3 (1 j1 j0 ) A2 (0 j1 j0 ) A2 (1 j1 j0 ).
195
表示为十进制,有
2
A
(
j
)
A
(
j
)
A
(
j
2
),
3
2
2
2
2
A3 ( j 2 ) A2 ( j ) A2 ( j 2 )
( j 0,1,2,3).
(6.15)
根据公式(6.13),(6.14),(6.15),由
A0 (k ) x(k ) xk
(k 0,1, ,7)
逐次计算到 A3 ( j ) c j ( j 0,1, ,7), 见表3-5(略).
从表3-5中看到计算全部8个 c j只用8次乘法运算和24次
加法运算.
196
上面推导的 N 23 的计算公式可类似地推广到 N 2 p
的情形.
根据公式(6.13),(6.14),(6.15),一般情况的FFT
计算公式如下:
q
q 1 j ) A
q 1 j 2 p 1 ),
q 1 (k 2
Aq (k 2 j ) Aq 1 (k 2
Aq (k 2q j 2q 1 )
(6.16)
q 1
[ Aq 1 (k 2 q 1 j ) Aq 1 ( k 2 q 1 j 2 p 1 )]wk 2 ,
其中 q 1,, p; k 0,1,,2 p q 1; j 0,1,,2q 1 1.
Aq 括号内的数代表它的位置,在计算机中代表存放数的地址.
197
一组 Aq占用 N个复数单元,计算时需给出两组单元,
从 A0 (m)( m 0,1,, N 1) 出发, q 由 1到 p算到
Ap ( j ) c j ( j 0,1, , N 1), 即为所求.
计算过程中只要按地址号存放 Aq , 则最后得到的 Ap ( j )
就是所求离散频谱的次序.
这个计算公式除了具有不倒地址的优点外,计算只有两
重循环.
p q
1,
外循环 q由 1计算到 p,内循环 k由 0计算到 2
j 由 0计算到 2 q 1 1, 更重要的是整个计算过程省计算量.
198
由公式看到算一个 Aq 共做 2 p q 2q /1 N / 2次复数乘法,
而最后一步计算 Ap 时,由于
k 2 p1
( N / 2 ) k (1) k (1) 0 1
(注意 q p 时 2 p q 1 0, 故 k 0 ),因此,总共要算
( p 1) N / 2 次复数乘法,它比直接用(6.9)需 N 2次乘法
快得多,计算量比值是 N : ( p 1) / 2.
当 N 210时比值是 1024 : 4.5 230 : 1, 它比一般FFT的
计算量( pN次乘法)也快一倍.
我们称(6.16)的计算公式为改进的FFT算法 .
199
例13
设 f ( x) x 4 3x 3 2 x 2 tan x( x 2) .给定
数据 {x j , f ( x j )}7j 0 , x j
解
j
确定三角插值多项式.
4
先将 [0,2] 变换为 [ , ] ,可令 y j ( x j 1) ,
故输入数据为 { y j , f j )}70 ,
f j f (1
yj
).
给定8个点可确定8个参数的4次三角插值多项式
3
1
S 4 ( y ) a0 (ak cos ky bk sin ky) a4 cos 4 y (6.17)
2
k 1
这里
2 7
2
a
f
cos
k 0,1, ,4,
k
j
8 kj,
8 j 0
7
2
b
f j sin 28 kj, k 1,2,3,
k
8 j 0
(6.18)
200
与(6.10)式比较可先计算
7
ck f j jk ,
j 0
这里 { f j }70代替(6.10)式的 {x j }70,
e
i 82
e
i 4
cos
4
isin
4
.
对每个 k 0,1,,4有
7
7
kj
1
1
1
1
i
ik ( 4 j )
k
ik
ik
4
ck ( 1) ck e
f je e
f je
4
4
4 j 0
4 j 0
1 7
j
j
f j (cosk ( ) isin k ( ))
4 j 0
4
4
1 7
f j (cosky j isin ky j ),
4 j 0
201
所以
(1) k
1
ak ibk
ck ck e ik ,
4
4
即
1
ak Re( ck e ik ),
4
1
bk Im( ck e ik ).
4
显然 b0 b4 0 ,用FFT算法求出 ck (k 0,1,2,3,4),也就
得到(6.18)式的系数,从而得到(6.17)式的4次三角
插值多项式
S 4 ( y ) 0.761979 0.771841cos y
0.386374 sin y 0.0173037 cos 2 y
0.0468750 sin 2 y 0.00686304 cos 3 y
0.0113738 sin 3 y 0.000578545 cos 4 y.
202
将 y ( x 1) 代入 S4 ( y) 可以得到 [0,2]上的三角多项式
S 4 ( x).
图3-6给出了 y f (x) 及 y S4 ( x) 的图形.
表3-6给出了在点 x j 0.125 j 0.25( j 0,1,,7) 处 f ( x j )
与 S 4 ( x j ) 的值.
203
图3-6
随机数据点上的二元拟合
204
205
206
207
208
209
210
211
212
213