上一页 下一页1
第二章 函数基本逼近 (一 )
-----插值逼近上一页 下一页2
当精确函数 y = f(x) 非常复杂或未知时,
在一系列节点
x0 … xn
处测得函数值
y0 = f(x0),…,yn = f(xn),
由此构造一个简单易算的近似函数 g(x)? f(x),
满足条件
g(xi) = f(xi) (i = 0,…,n)
这里的 g(x)称为 f(x) 的 插值函数 。
最常用的插值函数是 …?多项式上一页 下一页3
x0 x1 x2 x3 x4x
g(x)? f(x)
上一页 下一页4
§ 1 拉格朗日多项式 /* Lagrange Polynomial */
niyxP iin,...,0,)( ==
求 n 次多项式 使得 n
nn xaxaaxP=?10)(
条件,无重合节点,即 ji xx?ji?
n = 1 已知 x0,x1 ; y0,y1,求 xaaxP
101 )(?=
使得
111001 )(,)( yxPyxP ==
可见 P1(x) 是过 ( x0,y0 ) 和 ( x1,y1 ) 两点的直线。
)()( 0
01
01
01 xxxx
yyyxP -
-
-?=
10
1
xx
xx
-
-
01
0
xx
xx
-
-= y
0 + y1
l0(x) l1(x)
== 1
0
)(
i ii
yxl
称为 拉氏基函数 /* Lagrange Basis */,
满足条件 li(xj)=?ij /* Kronecker Delta */
上一页 下一页5
n? 1 希望找到 li(x),i = 0,…,n 使得 li(xj)=?ij ;然后令
=
=
n
i
iin yxlxP
0
)()(,则显然有 Pn(xj) = yj 。
li(x) 每个 li 有 n 个根 x0 … xi … xn
=
-=---= n
j
j? i
jiniii xxCxxxxxxCxl
0
0 )())...()...(()(
-== j? i
ji
iii xxCxl )( 11)(
=
-
-= n
j
ij ji
j
i xx
xxxl
0
)(
)()(
=
= n
i
iin yxlxL
0
)()(
Lagrange
Polynomial
与 有关,而与 无关节点 f
上一页 下一页6
定理 (唯一性 ) 满足 的 n 阶插值多项式是唯一存在的。
niyxP ii,...,0,)( ==
证明,( p 25 利用 Vandermonde 行列式 论证 )
反证:若不唯一,则除了 Ln(x) 外还有另一 n 阶多项式 Pn(x) 满足 Pn(xi) = yi 。
考察 则 Qn 的阶数,)()()( xLxPxQ
nnn -=
n
而 Qn 有 个不同的根n + 1 x0 … xn
注,若不将多项式次数限制为 n,则插值多项式 不唯一 。
例如 也是一个插值多项式,其中 可以是任意多项式。
=
-?= n
i
in xxxpxLxP
0
)()()()(
)(xp
上一页 下一页7
插值余项 /* Remainder */
设节点
)1(?nf 在 [a,b]内存在,考察截断误差
)()()( xLxfxR nn -=
],[ baCf n?bxxxa
n10
,且 f 满足条件,
Rolle’s Theorem,若 充分光滑,,则存在 使得 。
)(x? 0)()( 10 == xx
),( 10 xx 0)( =
推广,若 0)()()(
210 === xxx ),(),,( 211100 xxxx
使得 0)()(
10 =?= ),( 10
使得 0)( =
0)()( 0 === nxx
存在 ),( ba 使得 0)()( = n
Rn(x) 至少有 个根n+1?
=
-=
n
i
in xxxKxR
0
)()()(
任意固定 x? xi (i = 0,…,n),考察?
=
--=
n
i
ixtxKtRnt
0
)()()()(?
(t)有 n+2 个不同的根 x0 … xn x ),(,0)()1( ba
xxn?=
!)1()()()1(?-? nxKR xnn?
注意这里是对 t 求导
=?-- !)1)(()()( )1()1( nxKLf xnnxn
!)1(
)()( )1(
=
n
fxK xn?
=
-?=
n
i
i
x
n
n xxn
fxR
0
)1(
)(!)1( )()(?
上一页 下一页8
注,?通常不能确定?x,而是估计,?x?(a,b)
将 作为误差估计上限。
1)1( )( nn Mxf
=
-? n
i
in xxn
M
0
1 ||)!1(
当 f(x) 为任一个次数? n 的 多项式 时,,
可知,即插值多项式对于次数?n 的 多项式是 精确 的。
0)()1( xf n
0)(?xR n
Quiz:给定 xi = i +1,i = 0,1,2,3,4,5,下面哪个是 l2(x)的图像?
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
A B C
上一页 下一页9
例,已知
2 33s i n,214s i n,216s i n ===
分别利用 sin x 的 1次,2次 Lagrange 插值计算 sin 50?
并估计误差。
解:
0x 1x 2x
185500?=
n = 1 分别利用 x0,x1 以及 x1,x2 计算
4,6 10
== xx?利用
216/4/ 6/214/6/ 4/)(1?----= xxxL
这里 )
3,6(,s in)(,s in)( )2(-== xxxfxxf
而
)4)(6(!2 )()(,2 3s i n21 )2(1 --= xxfxR xx
0 0 7 6 2.0)185(0 1 3 1 9.0 1 --?R?sin 50?= 0.7660444…
)185(50sin 10L 0.77614
外推 /* extrapolation */ 的实际误差?-0.01010
3,4 21 == xx
利用 sin 50 0.76008,00660.0
185~00538.0 1R
内插 /* interpolation */ 的实际误差?0.00596
内插通常优于外推。选择要计算的 x 所在的区间的端点,插值效果较好。
上一页 下一页10
n = 2
2
3
))((
))((
2
1
))((
))((
21))((
))(()(
4363
46
3464
36
3646
342?-- ---- ---- --=
xxxxxxxL
)185(50sin 20L 0.76543
2 3c o s21;)3)(4)(6(!3
c o s)(
2---
-=
xx xxxxR
00 0 77.018500 0 44.0 2R?sin 50?= 0.7660444…
2次插值的实际误差?0.00061
高次插值通常优于低次插值但绝对不是次数越高就越好,……
上一页 下一页11
§ 2 牛顿插值 /* Newton’s Interpolation */
Lagrange 插值虽然易算,但若要增加一个节点时,
全部基函数 li(x) 都需重新算过。
将 Ln(x) 改写成,..))(()( 102010?--?-? xxxxaxxaa
)),,,(( 10 ---? nn xxxxa 的形式,希望每加一个节点时,
只附加一项 上去即可。
差商 (亦称均差 ) /* divided difference */
),()()(],[ ji
ji
ji
ji xxjixx
xfxfxxf
-
-=
1阶差商 /* the 1st
divided difference of f
w.r.t,xi and xj */
)(],[],[],,[ kixx xxfxxfxxxf
ki
kjji
kji?-
-= 2阶差商上一页 下一页12
1
11010
10
1110
10
],,...,[],,...,[
],,...,[],...,,[],...,[
--
-
-=
-
-=
kk
kkkk
k
kkk
k
xx
xxxfxxxf
xx
xxxfxxxfxxf
(k+1)阶差商:
=
= k
i ik
i
k x
xfxxf
0 1
0 )(
)(],.,,,[
事实上其中
,)()(
0
1?
=
-=
k
i
ik xxx
=
-=?
k
ij
j
jiik xxx
0
1 )()(?
差商的值与 xi 的顺序无关!
上一页 下一页13
牛顿插值 /* Newton’s Interpolation */
],[)()()( 000 xxfxxxfxf -?=
],,[)(],[],[ 101100 xxxfxxxxfxxf -?=
],...,,[)(],...,[],...,,[ 0010 nnnn xxxfxxxxfxxxf -?=-
)),,,((...))(()()( 10102010 -----?-?= nnn xxxxaxxxxaxxaaxN
1
2
… … … …
n-1
1 + (x - x0)? 2 + … … + ( x - x0)…( x - xn-1)? n-1
...))(](,,[)](,[)()( 102100100?--?-?= xxxxxxxfxxxxfxfxf
)),,,(](,...,[ 100 ---? nn xxxxxxf
))(),,,(](,...,,[ 100 nnn xxxxxxxxxf ---? -
Nn(x) R
n(x)
ai = f [ x0,…,xi ]
上一页 下一页14
注,?由 唯一性可知 Nn(x)? Ln(x),只是算法不同,故其余项也相同,即
)(!)1( )()(],...,,[ 1
)1(
10 xn
fxxxxf
k
x
n
kn?
=?
),(,! )(],...,[ m a xm i n)(0 xxkfxxf kk?=
实际计算过程为
f (x0)
f (x1)
f (x2)
…
f (xn-1)
f (xn)
f [x0,x1]
f [x1,x2]
… …
… …
f [xn-1,xn]
f [x0,x1,x2]
… …
… …
f [xn-2,xn-1,xn] f [x0,…,xn]
f (xn+1) f [xn,xn+1] f [xn-1,xn,xn+1] f [x1,…,xn+1] f [x0,…,xn+1]
上一页 下一页15
等距节点公式 /* Formulae with Equal Spacing */
向前差分
/* forward
difference */
iii fff -=1
i
k
i
k
i
k
i
k ffff 1
1
11 )( -
--?-?==?
向后差分
/* backward
difference */ 111 ---?-?=? ikikik fff
i-1ii fff -=?
中心差分
/* centered
difference */
2121
11 --?- -= ikikik fff
其中 )(
221 hii xff?=?
当节点 等距 分布时,),...,0(0 nihixx i =?=
上一页 下一页16
差分的重要性质:
线性:例如
gbfaxgbxfa= ))()((
若 f (x)是 m 次多项式,则 是 次多项式,而
)0()( mkxfk km-
)(0)( mkxfk?=?
差分值可由函数值算出:
=
-?-=?
n
j
jkn
j
k
n f
j
nf
0
)1(?
= -?
--=? n
j njk
jn
k
n f
j
nf
0
)1(
!
)1),,,(1(
j
jnnn
j
n?--=
其中 /* binomial coefficients */
函数值可由差分值算出,kj
n
j
kn fj
nf
=
=
0
k
k
k hk
fxxf
!],...,[
0
0
=
k
n
k
knnn hk
fxxxf
!],...,,[ 1
=
-- k
k
k
h
ff 0)( )(?=?
由 Rn 表达式上一页 下一页17
牛顿公式 )),,,(](,...,[...)](,[)()( 1000100 ----?= nnn xxxxxxfxxxxfxfxN
牛顿前差公式 /* Newton’s forward-difference formula */
牛顿后差公式 /* Newton’s backward-difference formula */
将节点顺序倒臵:
)),,,(](,...,[...)](,[)()( 101 xxxxxxfxxxxfxfxN nnnnnnn ---?= -
设 htxx?=
0
,则 )()()( 0
0
0 xfk
thtxNxN kn
k
nn?=?=?
=
),(,)),,,(1()!1( )()( 01)1( nnnn xxhntttnfxR?--?=
设 htxx
n?=
,则 )()1()()(
0
n
k
n
k
k
nnn xfk
thtxNxN?--=?=?
=
注,一般当 x 靠近 x0 时用前插,靠近 xn 时用后插,故两种公式亦称为 表初公式 和 表末公式 。
上一页 下一页18
§ 3 埃尔米特插值 /* Hermite Interpolation */
不仅要求函数值重合,而且要求若干阶 导数 也重合。
即:要求插值函数? (x) 满足?(xi) = f (xi),?’ (xi) = f ’ (xi),
…,?(mi) (xi) = f (mi) (xi).
注,?N 个条件可以确定 阶多项式。N - 1
要求在 1个节点 x0 处直到 m0阶导数都重合的插值多项式即为 Taylor多项式
0
0 )(
!
)(...))(()()(
00 0
)(
000
mm xxm xfxxxfxfx --=?
其余项为
)1(
0
0
)1(
0
0 )(
)!1(
)()()()( -
=-=
mm xx
m
fxxfxR
一般只考虑 f 与 f ’的值。
上一页 下一页19
例,设 x0? x1? x2,已知 f(x0),f(x1),f(x2) 和 f ’(x1),求多项式 P(x)
满足 P(xi) = f (xi),i = 0,1,2,且 P’(x1) = f ’(x1),并估计误差。
模仿 Lagrange 多项式的思想,设解,首先,P的阶数 = 3
= 2 13 )()()()()(
=0i ii
xhx1f ’xhxfxP?
h0(x) 有根 x1,x2,且 h0’(x1) = 0? x1 是重根。 )()()( 2
2
100 xxxxCxh --=
又,h0(x0) = 1? C0
)()(
)()()(
20
2
10
2
2
10
xxxx
xxxxxh
--
--=
h2(x)
h1(x) 有根 x0,x2? ))()(()( 201 xxxxBAxxh --?=
由余下条件 h1(x1) = 1 和 h1’(x1) = 0 可解。
与 h0(x) 完全类似。
(x)?h1 有根 x0,x1,x2h
1 ))()(()( 2101 xxxxxxCx ---=?
h1又,’(x1) = 1? C1 可解。
其中 hi(xj) =?ij,hi’(x1) = 0,(xi) = 0,’(x1) = 1?h1?h1
),())()(()()()( 221033 xxxxxxxKxPxfxR ---=-= !4 )()( )4(fxK?=
与 Lagrange 分析完全类似上一页 下一页20
一般地,已知 x0,…,xn 处有 y0,…,yn 和 y0’,…,yn’,求 H2n+1(x)
满足 H2n+1(xi) = yi,H’2n+1(xi) = yi’。
解,设= n i )()()(
=0i
i xhxhyixH2n+1
n
=0i
yi’
其中 hi(xj) =?ij,hi’(xj) = 0,(xj) = 0,’(xj) =?ij?hi?hi
hi(x) 有根 x0,…,xi,…,xn且都是 2重根 )()()( 2 xlBxAxh iiii?=
-
-=
ij ji
j
i xx
xxxl
)(
)()(
由余下条件 hi(xi) = 1 和 hi’(xi) = 0 可解 Ai 和 Bi?
)()])((21[)( 2 xlxxxlxh iiiii -?-=
(x)?hi 有根 x0,…,xn,除了 xi外都是 2重根hi )()( ii li2(x)xxCx -=
h
i又,’(xi) = 1? Ci = 1
h
i )(x )( i li2(x)xx-=
设 ],[,..,210 baCfbxxxa nn?== 则 2
0
)22(
)()!22( )()( -?=?
=
n
i
i
x
n
n xxn
fxR?
这样的 Hermite 插值唯一上一页 下一页21
Quiz,给定 xi = i +1,i = 0,1,2,3,4,5,下面哪个是 h2(x)的图像?
x0
-
-1
0.5
1 2 3 4 5 6
y
x
y
0
-
-
-1
0.5
1 2 3 4 5 6
斜率 =1?
求 Hermite多项式的基本步骤:
写出相应于条件的 hi(x),hi(x) 的组合形式;?
对每一个 hi(x),hi(x) 找出尽可能多的条件给出的根;?
根据多项式的总阶数和根的个数写出表达式;
根据尚未利用的条件解出表达式中的待定系数;
最后完整写出 H(x)。
上一页 下一页22
§ 4 分段低次插值 /* piecewise polynomial approximation */
例,在 [-5,5]上考察 的 Ln(x)。取
21
1)(
xxf?= ),.,,,0(105 niinx i =?-=
-5 -4 -3 -2 -1 0 1 2 3 4 5-0.5
0
0.5
1
1.5
2
2.5
n 越大,
端点附近抖动越大,称为
Runge 现象
Ln(x)? f (x)?
分段 低次 插值上一页 下一页23
分段线性插值 /* piecewise linear interpolation */
在每个区间 上,用 1阶多项式 (直线 ) 逼近 f (x):],[ 1?ii xx
11111 )()( -
-?
-
-=?
iii iiii i yxx
xxy
xx
xxxPxf ],[f o r
1 ii xxx
记,易证:当 时,||m a x
1 ii xxh -=? 0?h )()(1 xfxP h?
一致失去了原函数的光滑性。
上一页 下一页24
分段 Hermite插值 /* Hermite piecewise polynomials */
给定
nnn yyyyxx,...,;,...,;,...,000
在 上利用两点的 y 及 y’ 构造 3次 Hermite函数],[ 1?ii xx
导数一般不易得到。
上一页 下一页25
§ 5 三次样条 /* Cubic Spline */
定义 设 。 三次样条函数,
且在每个 上为 三次多项式 /* cubic polynomial */。 若它同时还满足,则称为 f 的 三次样条插值函数 /* cubic spline interpolant */.
bxxxa n ==,..10 ],[)( 2 baCxS?
],[ 1?ii xx
),...,0(),()( nixfxS ii ==
注,三次样条与分段 Hermite 插值的根本区别在于 S(x)自身光滑,不需要知道 f 的导数值(除了在 2个端点可能需要);而 Hermite插值依赖于 f 在所有插值点的导数值。
f(x) H(x)
S(x)
上一页 下一页26
构造三次样条插值函数的 三弯矩法
/* method of bending moment */
在 上,记],[
1 jj xx -,1--= jjj xxh ],[for )()( 1
][
jj
j xxxxSxS
-?=
对每个 j,此为 3次多项式则 S[j]”(x) 为 次多项式,需 个点的值确定之。1 2
设 S[j]”(xj-1) = Mj-1,S[j]”(xj) = Mj
对应力学中的 梁弯矩,故名对于 x?[xj-1,xj ] 可 得到
S[j]”(x) =
j
j
j
j
j
j h
xxM
h
xxM 1
1
-
-
-?-
积分 2次,可得 S[j]’(x) 和 S[j](x),
j
j
j
j
j
j
j Ah
xxM
h
xxM?-?-- -
-- 2
)(
2
)( 21
1
2
1S
[j]’(x) =
jj
j
j
j
j
j
j BxAh
xxM
h
xxM-?- -
- 6
)(
6
)( 313
1S
[j](x) =
利用已知
S[j](xj-1) = yj-1
S[j](xj) = yj
可解上一页 下一页27
j
jj
j
jj
j h
MM
h
yyA
6
11 -- ---=
j
j
j
j
j
j
j
j
j
jjj h
xxhMy
h
xxhMyBxA 1221
1 )6()6(
--
-
--?--=?
下面解决 Mj,利用 S’在 xj 的 连续性
[xj-1,xj ],S[j]’(x) = jjjjj
j
j
j
j
j
j h
MMxxf
h
xxM
h
xxM
6],[2
)(
2
)( 1
1
2
1
2
1
-
-
-
-
--?-?--
1
1
1
1
2
1
1
2
1
6],[2
)(
2
)(
--?-?--
j
jj
jj
j
j
j
j
j
j h
MMxxf
h
xxM
h
xxM[x
j,xj+1],S[j+1]’(x) =
利用 S[j]’(xj) = S[j+1]’(xj),合并关于 Mj-1,Mj,Mj+1的同类项,并记,,,整理后得到:
1
1
jj
j
j hh
h
=l 1 jj -= lm ]),[],[(
6
11
1
jjjj
jj
j xxfxxfhhg -?
-?=
2 11 gMMM jjjjjj =- lm
j?1 n-1
即:有 个未知数,个方程。n-1n+1
=
--- 1
1
0
11
11
2
2
n
n
nn
g
g
M
M
lm
lm
还需 2个 边界条件 /* boundary conditions */
上一页 下一页28
第 1类边条件 /* clamped boundary */,S’(a) = y0’,S’(b) = yn’
[a,x1 ],S[1]’(x) = 10110
1
2
1
1
2
1
0 6],[2
)(
2
)( hMMxxf
h
axM
h
xxM --?-?--
010
1
10 )],[(
62 gy
0’xxfhMM =-=?
nnn
n
nn gxxfyn’hMM =-=? -- ]),[(
62
11
类似地利用 [ xn-1,b ]
上的 S[n]’(x)
第 2类边条件,S”(a) = y0” = M0,S”(b) = yn” = Mn
这时:
nnn ygyg==== 2,0;2,0 000 ml
特别地,M0 = Mn = 0 称为 自由边界 /* free boundary */,对应的样条函数称为 自然样条 /* Natural Spline */。
第 3类边条件 /* periodic boundary */,
当 f 为 周期函数 时,
yn = y0,S’(a+) = S’(b-)
M0 = Mn
=
--
nnnn
nn
g
g
M
M 11
11
22
11
2
2
2
2
ml
lm
lm
ml
上一页 下一页29
注,?另有 三转角法 得到样条函数,即设 S[j]’(xj) = mj,则易知 [xj-1,xj ]上的 S[j](x) 就是 Hermite函数 。再利用 S”
的连续性,可导出关于 mj 的方程组,加上边界条件即可解。
Cubic Spline 由 boundary conditions 唯一 确定。
收敛性,若,且,则],[ baCf C
hh iim i nm a x一致
S(x)?f(x) 0m a x?
ihas
即,提高精度只须 增加节点,而无须提高样条阶数。
稳定性,只要边条件保证 | m0 |,| l0 |,| mn |,| l n | < 2,
则方程组系数阵为 SDD阵,保证数值稳定。
上一页 下一页30
Sketch of the Algorithm,Cubic Spline
① 计算 mj,l j,gj ; ② 计算 Mj (追赶法等 ) ;
③ 找到 x 所在区间 ( 即找到相应的 j ) ;
④ 由该区间上的 S[j](x) 算出 f(x) 的近似值。
插值法小结
Lagrange,给出 y0 … yn,选基函数 li(x),其次数为节点数 –1。
Newton? Ln(x),只是形式不同;节点等距或渐增节点时方便处理。
Hermite,给出 yi 及 yi ’,选 hi(x) 及 hi(x) 。
Spline:分段低次,自身光滑,f 的导数只在边界给出。
上一页 下一页31
作业,p.56-59
1题 (1),2题,4题,
5题,7题,9题,
10题,13题,20题。
上一页 下一页32
当精确函数 y = f(x) 非常复杂或未知时,在一系列节点 x0 … xn 处测得函数值 y0 = f(x0),…
yn = f(xn),由此构造一个简单易算的近似函数 g(x)? f(x),满足条件 g(xi) = f(xi) (i = 0,…
n)。这里的 g(x)称为 f(x) 的 插值函数 。最常用的插值函数是 …?多项式
x0 x1 x2 x3 x4x
g(x)? f(x)
上一页 下一页33
§ 1 拉格朗日多项式 /* Lagrange Polynomial */
niyxP iin,...,0,)( ==
求 n 次多项式 使得 n
nn xaxaaxP=?10)(
条件,无重合节点,即 ji xx?ji?
n = 1 已知 x0,x1 ; y0,y1,求 xaaxP
101 )(?=
使得
111001 )(,)( yxPyxP ==
可见 P1(x) 是过 ( x0,y0 ) 和 ( x1,y1 ) 两点的直线。
)()( 0
01
01
01 xxxx
yyyxP -
-
-?=
10
1
xx
xx
-
-
01
0
xx
xx
-
-= y
0 + y1
l0(x) l1(x)
== 1
0
)(
i ii
yxl
称为 拉氏基函数 /* Lagrange Basis */,
满足条件 li(xj)=?ij /* Kronecker Delta */
上一页 下一页34
n? 1 希望找到 li(x),i = 0,…,n 使得 li(xj)=?ij ;然后令
=
=
n
i
iin yxlxP
0
)()(,则显然有 Pn(xi) = yi 。
li(x) 每个 li 有 n 个根 x0 … xi … xn
=
-=---= n
j
j? i
jiniii xxCxxxxxxCxl
0
0 )())...()...(()(
-== j? i
ji
iii xxCxl )( 11)(
=
-
-= n
j
ij ji
j
i xx
xxxl
0
)(
)()(
=
= n
i
iin yxlxL
0
)()(
Lagrange
Polynomial
与 有关,而与 无关节点 f
上一页 下一页35
定理 (唯一性 ) 满足 的 n 阶插值多项式是唯一存在的。
niyxP ii,...,0,)( ==
证明,( p 25 利用 Vandermonde 行列式 论证 )
反证:若不唯一,则除了 Ln(x) 外还有另一 n 阶多项式 Pn(x) 满足 Pn(xi) = yi 。
考察 则 Qn 的阶数,)()()( xLxPxQ
nnn -=
n
而 Qn 有 个不同的根n + 1 x0 … xn
注,若不将多项式次数限制为 n,则插值多项式 不唯一 。
例如 也是一个插值多项式,其中 可以是任意多项式。
=
-?= n
i
in xxxpxLxP
0
)()()()(
)(xp
上一页 下一页36
插值余项 /* Remainder */
设节点
)1(?nf 在 [a,b]内存在,考察截断误差
)()()( xLxfxR nn -=
],[ baCf n?bxxxa
n10
,且 f 满足条件,
Rolle’s Theorem,若 充分光滑,,则存在 使得 。
)(x? 0)()( 10 == xx
),( 10 xx 0)( =
推广,若 0)()()(
210 === xxx ),(),,( 211100 xxxx
使得 0)()(
10 =?= ),( 10
使得 0)( =
0)()( 0 === nxx
存在 ),( ba 使得 0)()( = n
Rn(x) 至少有 个根n+1?
=
-=
n
i
in xxxKxR
0
)()()(
任意固定 x? xi (i = 0,…,n),考察?
=
--=
n
i
ixtxKtRnt
0
)()()()(?
(x)有 n+2 个不同的根 x0… xn x ),(,0)()1( ba
xxn?=
!)1()()()1(?-? nxKR xnn?
注意这里是对 t 求导
=?-- !)1)(()()( )1()1( nxKLf xnnxn
!)1(
)()( )1(
=
n
fxK xn?
=
-?=
n
i
i
x
n
n xxn
fxR
0
)1(
)(!)1( )()(?
上一页 下一页37
注,?通常不能确定?x,而是估计,?x?(a,b)
将 作为误差估计上限。
1)1( )( nn Mxf
=
-? n
i
in xxn
M
0
1 ||)!1(
当 f(x) 为任一个次数? n 的 多项式 时,,
可知,即插值多项式对于次数?n 的 多项式是 精确 的。
0)()1( xf n
0)(?xR n
Quiz:给定 xi = i +1,i = 0,1,2,3,4,5,下面哪个是 l2(x)的图像?
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
A B C
上一页 下一页38
例,已知
2 33s i n,214s i n,216s i n ===
分别利用 sin x 的 1次,2次 Lagrange 插值计算 sin 50?
并估计误差。
解:
0x 1x 2x
185500?=
n = 1 分别利用 x0,x1 以及 x1,x2 计算
4,6 10
== xx?利用
216/4/ 6/214/6/ 4/)(1?----= xxxL
这里 )
3,6(,s in)(,s in)( )2(-== xxxfxxf
而
)4)(6(!2 )()(,2 3s i n21 )2(1 --= xxfxR xx
0 0 7 6 2.0)185(0 1 3 1 9.0 1 --?R?sin 50?= 0.7660444…
)185(50sin 10L 0.77614
外推 /* extrapolation */ 的实际误差?-0.01010
3,4 21 == xx
利用 sin 50 0.76008,00660.0
185~00538.0 1R
内插 /* interpolation */ 的实际误差?0.00596
内插通常优于外推。选择要计算的 x 所在的区间的端点,插值效果较好。
上一页 下一页39
n = 2
2
3
))((
))((
2
1
))((
))((
21))((
))(()(
4363
46
3464
36
3646
342?-- ---- ---- --=
xxxxxxxL
)185(50sin 20L 0.76543
2 3c o s21;)3)(4)(6(!3
c o s)(
2---
-=
xx xxxxR
00 0 77.018500 0 44.0 2R?sin 50?= 0.7660444…
2次插值的实际误差?0.00061
高次插值通常优于低次插值但绝对不是次数越高就越好,嘿嘿 ……
上一页 下一页40
§ 2 牛顿插值 /* Newton’s Interpolation */
Lagrange 插值虽然易算,但若要增加一个节点时,
全部基函数 li(x) 都需重新算过。
将 Ln(x) 改写成,..))(()( 102010?--?-? xxxxaxxaa
)),,,(( 10 ---? nn xxxxa 的形式,希望每加一个节点时,
只附加一项 上去即可。
差商 (亦称均差 ) /* divided difference */
),()()(],[ ji
ji
ji
ji xxjixx
xfxfxxf
-
-=
1阶差商 /* the 1st
divided difference of f
w.r.t,xi and xj */
)(],[],[],,[ kixx xxfxxfxxxf
ki
kjji
kji?-
-= 2阶差商上一页 下一页41
1
11010
10
1110
10
],,...,[],,...,[
],,...,[],...,,[],...,[
--
-
-=
-
-=
kk
kkkk
k
kkk
k
xx
xxxfxxxf
xx
xxxfxxxfxxf
(k+1)阶差商:
=
= k
i ik
i
k x
xfxxf
0 1
0 )(
)(],.,,,[
事实上其中
,)()(
0
1?
=
-=
k
i
ik xxx
=
-=?
k
ij
j
jiik xxx
0
1 )()(?
差商的值与 xi 的顺序无关!
上一页 下一页42
牛顿插值 /* Newton’s Interpolation */
],[)()()( 000 xxfxxxfxf -?=
],,[)(],[],[ 101100 xxxfxxxxfxxf -?=
],...,,[)(],...,[],...,,[ 0010 nnnn xxxfxxxxfxxxf -?=-
)),,,((...))(()()( 10102010 -----?-?= nnn xxxxaxxxxaxxaaxN
1
2
… … … …
n-1
1 + (x - x0)? 2 + … … + ( x - x0)…( x - xn-1)? n-1
...))(](,,[)](,[)()( 102100100?--?-?= xxxxxxxfxxxxfxfxf
)),,,(](,...,[ 100 ---? nn xxxxxxf
))(),,,(](,...,,[ 100 nnn xxxxxxxxxf ---? -
Nn(x) R
n(x)
ai = f [ x0,…,xi ]
上一页 下一页43
注,?由 唯一性可知 Nn(x)? Ln(x),只是算法不同,故其余项也相同,即
)(!)1( )()(],...,,[ 1
)1(
10 xn
fxxxxf
k
x
n
kn?
=?
),(,! )(],...,[ m a xm i n)(0 xxkfxxf kk?=
实际计算过程为
f (x0)
f (x1)
f (x2)
…
f (xn-1)
f (xn)
f [x0,x1]
f [x1,x2]
… …
… …
f [xn-1,xn]
f [x0,x1,x2]
… …
… …
f [xn-2,xn-1,xn] f [x0,…,xn]
f (xn+1) f [xn,xn+1] f [xn-1,xn,xn+1] f [x1,…,xn+1] f [x0,…,xn+1]
上一页 下一页44
等距节点公式 /* Formulae with Equal Spacing */
向前差分
/* forward
difference */
iii fff -=1
i
k
i
k
i
k
i
k ffff 1
1
11 )( -
--?-?==?
向后差分
/* backward
difference */ 111 ---?-?=? ikikik fff
i-1ii fff -=?
中心差分
/* centered
difference */
2121
11 --?- -= ikikik fff
其中 )(
221 hii xff?=?
当节点 等距 分布时,),...,0(0 nihixx i =?=
上一页 下一页45
差分的重要性质:
线性:例如
gbfaxgbxfa= ))()((
若 f (x)是 m 次多项式,则 是 次多项式,而
)0()( mkxfk km-
)(0)( mkxfk?=?
差分值可由函数值算出:
=
-?-=?
n
j
jkn
j
k
n f
j
nf
0
)1(?
= -?
--=? n
j njk
jn
k
n f
j
nf
0
)1(
!
)1),,,(1(
j
jnnn
j
n?--=
其中 /* binomial coefficients */
函数值可由差分值算出,kj
n
j
kn fj
nf
=
=
0
k
k
k hk
fxxf
!],...,[
0
0
=
k
n
k
knnn hk
fxxxf
!],...,,[ 1
=
-- k
k
k
h
ff 0)( )(?=?
由 Rn 表达式上一页 下一页46
牛顿公式 )),,,(](,...,[...)](,[)()( 1000100 ----?= nnn xxxxxxfxxxxfxfxN
牛顿前差公式 /* Newton’s forward-difference formula */
牛顿后差公式 /* Newton’s backward-difference formula */
将节点顺序倒臵:
)),,,(](,...,[...)](,[)()( 101 xxxxxxfxxxxfxfxN nnnnnnn ---?= -
设 htxx?=
0
,则 )()()( 0
0
0 xfk
thtxNxN kn
k
nn?=?=?
=
),(,)),,,(1()!1( )()( 01)1( nnnn xxhntttnfxR?--?=
设 htxx
n?=
,则 )()1()()(
0
n
k
n
k
k
nnn xfk
thtxNxN?--=?=?
=
注,一般当 x 靠近 x0 时用前插,靠近 xn 时用后插,故两种公式亦称为 表初公式 和 表末公式 。
上一页 下一页47
§ 3 埃尔米特插值 /* Hermite Interpolation */
不仅要求函数值重合,而且要求若干阶 导数 也重合。
即:要求插值函数? (x) 满足?(xi) = f (xi),?’ (xi) = f ’ (xi),
…,?(mi) (xi) = f (mi) (xi).
注,?N 个条件可以确定 阶多项式。N - 1
要求在 1个节点 x0 处直到 m0阶导数都重合的插值多项式即为 Taylor多项式
0
0 )(
!
)(...))(()()(
00 0
)(
000
mm xxm xfxxxfxfx --=?
其余项为
)1(
0
0
)1(
0
0 )(
)!1(
)()()()( -
=-=
mm xx
m
fxxfxR
一般只考虑 f 与 f ’的值。
上一页 下一页48
例,设 x0? x1? x2,已知 f(x0),f(x1),f(x2) 和 f ’(x1),求多项式 P(x)
满足 P(xi) = f (xi),i = 0,1,2,且 P’(x1) = f ’(x1),并估计误差。
模仿 Lagrange 多项式的思想,设解,首先,P的阶数 = 3
= 2 13 )()()()()(
=0i ii
xhx1f ’xhxfxP?
h0(x) 有根 x1,x2,且 h0’(x1) = 0? x1 是重根。 )()()( 2
2
100 xxxxCxh --=
又,h0(x0) = 1? C0
)()(
)()()(
20
2
10
2
2
10
xxxx
xxxxxh
--
--=
h2(x)
h1(x) 有根 x0,x2? ))()(()( 201 xxxxBAxxh --?=
由余下条件 h1(x1) = 1 和 h1’(x1) = 0 可解。
与 h0(x) 完全类似。
(x)?h1 有根 x0,x1,x2h
1 ))()(()( 2101 xxxxxxCx ---=?
h1又,’(x1) = 1? C1 可解。
其中 hi(xj) =?ij,hi’(x1) = 0,(xi) = 0,’(x1) = 1?h1?h1
),())()(()()()( 221033 xxxxxxxKxPxfxR ---=-= !4 )()( )4(fxK?=
与 Lagrange 分析完全类似上一页 下一页49
一般地,已知 x0,…,xn 处有 y0,…,yn 和 y0’,…,yn’,求 H2n+1(x)
满足 H2n+1(xi) = yi,H’2n+1(xi) = yi’。
解,设= n i )()()(
=0i
i xhxhyixH2n+1
n
=0i
yi’
其中 hi(xj) =?ij,hi’(xj) = 0,(xj) = 0,’(xj) =?ij?hi?hi
hi(x) 有根 x0,…,xi,…,xn且都是 2重根 )()()( 2 xlBxAxh iiii?=
-
-=
ij ji
j
i xx
xxxl
)(
)()(
由余下条件 hi(xi) = 1 和 hi’(xi) = 0 可解 Ai 和 Bi?
)()])((21[)( 2 xlxxxlxh iiiii -?-=
(x)?hi 有根 x0,…,xn,除了 xi外都是 2重根hi )()( ii li2(x)xxCx -=
h
i又,’(xi) = 1? Ci = 1
h
i )(x )( i li2(x)xx-=
设 ],[,..,210 baCfbxxxa nn?== 则 2
0
)22(
)()!22( )()( -?=?
=
n
i
i
x
n
n xxn
fxR?
这样的 Hermite 插值唯一上一页 下一页50
Quiz,给定 xi = i +1,i = 0,1,2,3,4,5,下面哪个是 h2(x)的图像?
x0
-
-1
0.5
1 2 3 4 5 6
y
x
y
0
-
-
-1
0.5
1 2 3 4 5 6
斜率 =1?
求 Hermite多项式的基本步骤:
写出相应于条件的 hi(x),hi(x) 的组合形式;?
对每一个 hi(x),hi(x) 找出尽可能多的条件给出的根;?
根据多项式的总阶数和根的个数写出表达式;
根据尚未利用的条件解出表达式中的待定系数;
最后完整写出 H(x)。
上一页 下一页51
§ 4 分段低次插值 /* piecewise polynomial approximation */
例,在 [-5,5]上考察 的 Ln(x)。取
21
1)(
xxf?= ),.,,,0(105 niinx i =?-=
-5 -4 -3 -2 -1 0 1 2 3 4 5-0.5
0
0.5
1
1.5
2
2.5
n 越大,
端点附近抖动越大,称为
Runge 现象
Ln(x)? f (x)?
分段 低次 插值上一页 下一页52
分段线性插值 /* piecewise linear interpolation */
在每个区间 上,用 1阶多项式 (直线 ) 逼近 f (x):],[ 1?ii xx
11111 )()( -
-?
-
-=?
iii iiii i yxx
xxy
xx
xxxPxf ],[f o r
1 ii xxx
记,易证:当 时,||m a x
1 ii xxh -=? 0?h )()(1 xfxP h?
一致失去了原函数的光滑性。
上一页 下一页53
分段 Hermite插值 /* Hermite piecewise polynomials */
给定
nnn yyyyxx,...,;,...,;,...,000
在 上利用两点的 y 及 y’ 构造 3次 Hermite函数],[ 1?ii xx
导数一般不易得到。
上一页 下一页54
§ 5 三次样条 /* Cubic Spline */
定义 设 。 三次样条函数,
且在每个 上为 三次多项式 /* cubic polynomial */。 若它同时还满足,则称为 f 的 三次样条插值函数 /* cubic spline interpolant */.
bxxxa n ==,..10 ],[)( 2 baCxS?
],[ 1?ii xx
),...,0(),()( nixfxS ii ==
注,三次样条与分段 Hermite 插值的根本区别在于 S(x)自身光滑,不需要知道 f 的导数值(除了在 2个端点可能需要);而 Hermite插值依赖于 f 在所有插值点的导数值。
f(x) H(x)
S(x)
上一页 下一页55
构造三次样条插值函数的 三弯矩法
/* method of bending moment */
在 上,记],[
1 jj xx -,1--= jjj xxh ],[for )()( 1
][
jj
j xxxxSxS
-?=
对每个 j,此为 3次多项式则 S[j]”(x) 为 次多项式,需 个点的值确定之。1 2
设 S[j]”(xj-1) = Mj-1,S[j]”(xj) = Mj
对应力学中的 梁弯矩,故名对于 x?[xj-1,xj ] 可 得到
S[j]”(x) =
j
j
j
j
j
j h
xxM
h
xxM 1
1
-
-
-?-
积分 2次,可得 S[j]’(x) 和 S[j](x),
j
j
j
j
j
j
j Ah
xxM
h
xxM?-?-- -
-- 2
)(
2
)( 21
1
2
1S
[j]’(x) =
jj
j
j
j
j
j
j BxAh
xxM
h
xxM-?- -
- 6
)(
6
)( 313
1S
[j](x) =
利用已知
S[j](xj-1) = yj-1
S[j](xj) = yj
可解上一页 下一页56
j
jj
j
jj
j h
MM
h
yyA
6
11 -- ---=
j
j
j
j
j
j
j
j
j
jjj h
xxhMy
h
xxhMyBxA 1221
1 )6()6(
--
-
--?--=?
下面解决 Mj,利用 S’在 xj 的 连续性
[xj-1,xj ],S[j]’(x) = jjjjj
j
j
j
j
j
j h
MMxxf
h
xxM
h
xxM
6],[2
)(
2
)( 1
1
2
1
2
1
-
-
-
-
--?-?--
1
1
1
1
2
1
1
2
1
6],[2
)(
2
)(
--?-?--
j
jj
jj
j
j
j
j
j
j h
MMxxf
h
xxM
h
xxM[x
j,xj+1],S[j+1]’(x) =
利用 S[j]’(xj) = S[j+1]’(xj),合并关于 Mj-1,Mj,Mj+1的同类项,并记,,,整理后得到:
1
1
jj
j
j hh
h
=l 1 jj -= lm ]),[],[(
6
11
1
jjjj
jj
j xxfxxfhhg -?
-?=
2 11 gMMM jjjjjj =- lm
j?1 n-1
即:有 个未知数,个方程。n-1n+1
=
--- 1
1
0
11
11
2
2
n
n
nn
g
g
M
M
lm
lm
还需 2个 边界条件 /* boundary conditions */
上一页 下一页57
第 1类边条件 /* clamped boundary */,S’(a) = y0’,S’(b) = yn’
[a,x1 ],S[1]’(x) = 10110
1
2
1
1
2
1
0 6],[2
)(
2
)( hMMxxf
h
axM
h
xxM --?-?--
010
1
10 )],[(
62 gy
0’xxfhMM =-=?
nnn
n
nn gxxfyn’hMM =-=? -- ]),[(
62
11
类似地利用 [ xn-1,b ]
上的 S[n]’(x)
第 2类边条件,S”(a) = y0” = M0,S”(b) = yn” = Mn
这时:
nnn ygyg==== 2,0;2,0 000 ml
特别地,M0 = Mn = 0 称为 自由边界 /* free boundary */,对应的样条函数称为 自然样条 /* Natural Spline */。
第 3类边条件 /* periodic boundary */,
当 f 为 周期函数 时,
yn = y0,S’(a+) = S’(b-)
M0 = Mn
=
--
nnnn
nn
g
g
M
M 11
11
22
11
2
2
2
2
ml
lm
lm
ml
上一页 下一页58
注,?另有 三转角法 得到样条函数,即设 S[j]’(xj) = mj,则易知 [xj-1,xj ]上的 S[j](x) 就是 Hermite函数 。再利用 S”
的连续性,可导出关于 mj 的方程组,加上边界条件即可解。
Cubic Spline 由 boundary conditions 唯一 确定。
收敛性,若,且,则],[ baCf C
hh iim i nm a x一致
S(x)?f(x) 0m a x?
ihas
即,提高精度只须 增加节点,而无须提高样条阶数。
稳定性,只要边条件保证 | m0 |,| l0 |,| mn |,| l n | < 2,
则方程组系数阵为 SDD阵,保证数值稳定。
上一页 下一页59
Sketch of the Algorithm,Cubic Spline
① 计算 mj,l j,gj ; ② 计算 Mj (追赶法等 ) ;
③ 找到 x 所在区间 ( 即找到相应的 j ) ;
④ 由该区间上的 S[j](x) 算出 f(x) 的近似值。
插值法小结
Lagrange,给出 y0 … yn,选基函数 li(x),其次数为节点数 –1。
Newton? Ln(x),只是形式不同;节点等距或渐增节点时方便处理。
Hermite,给出 yi 及 yi ’,选 hi(x) 及 hi(x) 。
Spline:分段低次,自身光滑,f 的导数只在边界给出。
上一页 下一页60
§ 5 Cubic Spline
Lab 1,Cubic Spline
Construct the cubic spline interpolant S for the function f,defined at
points,satisfying some given boundary conditions,Partition
an given interval into m equal-length subintervals,and approximate the
function values at the endpoints of these subintervals.
Input
There are several sets of inputs,For each set:
The 1st line contains an integer 20? n? 0 which is the number of interpolating
points - 1,n = -1 signals the end of file.
The 2nd line contains n+1 distinct real numbers,
The 3rd line contains n+1 real numbers,
The 4th line contains an integer Type and three real numbers s0,sn,and Fmax.
If Type = 1,the clamped boundary condition is given,and,
If Type = 2,the natural boundary condition is given,and,
No other values of Type will be given.
Fmax is the default value that S (x) will assume if x is out of the range [x0,xn].
The last line of a test case consists of two real numbers t0 and tm,and an integer m?
1,Here the interval is partitioned into m equal-lengh subintervals,
…,The numbers are separated by spaces and new lines.
nxxx,...,,10
)(,...),(),( 10 nxfxfxf
nxxx,..10
00 )( sxS =? nn sxS =? )(
00 )( sxS = nn sxS = )(
],[ 0 mtt ],[ 10 tt
],[ 1 mm tt -
上一页 下一页61
§ 5 Cubic Spline
Output (? represents a space)
For each test case,you are supposed to output the following information:
1,The set of coefficients of S(x) in the format:
where
for,
Each coefficient is to be printed as in C printf:
fprintf(outfile,"%12.8e? ",coefficient );
2,m+1 approximating values of f at in the format:
fprintf(outfile,"f(%12.8e)? =? %12.8e\n",t,S );
The outputs of two test cases must be seperated by a blank line.
nnnn dcba
dcba
dcba
..
.
..
.
..
.
..
.
2222
1111
31211][ )()()()()( --- -?-?-?== jjjjjjjj xxdxxcxxbaxSxS
],[ 1 jj xxx -?
mtt,..,,0
Sample Input
2
0.0? 1.0? 2.0
0.0? 1.0? 2.0
1? 1.0? 1.0? 0.0
0.0? 3.0? 2
2
–0.5? –0.25? 0.0
–0.0247500? 0.3349375? 1.1010000
2? 0.0? 0.0? 0.0
–1.0? 0.0? 4
–1
上一页 下一页62
§ 5 Cubic Spline
Sample Output
0.00000000e+000? 1.00000000e+000? 0.00000000e+000? 0.00000000e+000?
1.00000000e+000? 1.00000000e+000? 0.00000000e+000? 0.00000000e+000?
f(0.00000000e+000)? =? 0.00000000e+000
f(1.50000000e+000)? =? 1.50000000e+000
f(3.00000000e+000)? =? 0.00000000e+000
–2.47500000e–002? 1.03237500e+000? 0.00000000e+000? 6.50200000e+000?
3.34937500e–001? 2.25150000e+000? 4.87650000e+000? –6.50200000e+000?
f(–1.00000000e+000)? =? 0.00000000e+000
f(–7.50000000e–001)? =? 0.00000000e+000
f(–5.00000000e–001)? =? –2.47500000e–002
f(–2.50000000e–001)? =? 3.34937500e–001
f(0.00000000e+000)? =? 1.10100000e+000
第二章 函数基本逼近 (一 )
-----插值逼近上一页 下一页2
当精确函数 y = f(x) 非常复杂或未知时,
在一系列节点
x0 … xn
处测得函数值
y0 = f(x0),…,yn = f(xn),
由此构造一个简单易算的近似函数 g(x)? f(x),
满足条件
g(xi) = f(xi) (i = 0,…,n)
这里的 g(x)称为 f(x) 的 插值函数 。
最常用的插值函数是 …?多项式上一页 下一页3
x0 x1 x2 x3 x4x
g(x)? f(x)
上一页 下一页4
§ 1 拉格朗日多项式 /* Lagrange Polynomial */
niyxP iin,...,0,)( ==
求 n 次多项式 使得 n
nn xaxaaxP=?10)(
条件,无重合节点,即 ji xx?ji?
n = 1 已知 x0,x1 ; y0,y1,求 xaaxP
101 )(?=
使得
111001 )(,)( yxPyxP ==
可见 P1(x) 是过 ( x0,y0 ) 和 ( x1,y1 ) 两点的直线。
)()( 0
01
01
01 xxxx
yyyxP -
-
-?=
10
1
xx
xx
-
-
01
0
xx
xx
-
-= y
0 + y1
l0(x) l1(x)
== 1
0
)(
i ii
yxl
称为 拉氏基函数 /* Lagrange Basis */,
满足条件 li(xj)=?ij /* Kronecker Delta */
上一页 下一页5
n? 1 希望找到 li(x),i = 0,…,n 使得 li(xj)=?ij ;然后令
=
=
n
i
iin yxlxP
0
)()(,则显然有 Pn(xj) = yj 。
li(x) 每个 li 有 n 个根 x0 … xi … xn
=
-=---= n
j
j? i
jiniii xxCxxxxxxCxl
0
0 )())...()...(()(
-== j? i
ji
iii xxCxl )( 11)(
=
-
-= n
j
ij ji
j
i xx
xxxl
0
)(
)()(
=
= n
i
iin yxlxL
0
)()(
Lagrange
Polynomial
与 有关,而与 无关节点 f
上一页 下一页6
定理 (唯一性 ) 满足 的 n 阶插值多项式是唯一存在的。
niyxP ii,...,0,)( ==
证明,( p 25 利用 Vandermonde 行列式 论证 )
反证:若不唯一,则除了 Ln(x) 外还有另一 n 阶多项式 Pn(x) 满足 Pn(xi) = yi 。
考察 则 Qn 的阶数,)()()( xLxPxQ
nnn -=
n
而 Qn 有 个不同的根n + 1 x0 … xn
注,若不将多项式次数限制为 n,则插值多项式 不唯一 。
例如 也是一个插值多项式,其中 可以是任意多项式。
=
-?= n
i
in xxxpxLxP
0
)()()()(
)(xp
上一页 下一页7
插值余项 /* Remainder */
设节点
)1(?nf 在 [a,b]内存在,考察截断误差
)()()( xLxfxR nn -=
],[ baCf n?bxxxa
n10
,且 f 满足条件,
Rolle’s Theorem,若 充分光滑,,则存在 使得 。
)(x? 0)()( 10 == xx
),( 10 xx 0)( =
推广,若 0)()()(
210 === xxx ),(),,( 211100 xxxx
使得 0)()(
10 =?= ),( 10
使得 0)( =
0)()( 0 === nxx
存在 ),( ba 使得 0)()( = n
Rn(x) 至少有 个根n+1?
=
-=
n
i
in xxxKxR
0
)()()(
任意固定 x? xi (i = 0,…,n),考察?
=
--=
n
i
ixtxKtRnt
0
)()()()(?
(t)有 n+2 个不同的根 x0 … xn x ),(,0)()1( ba
xxn?=
!)1()()()1(?-? nxKR xnn?
注意这里是对 t 求导
=?-- !)1)(()()( )1()1( nxKLf xnnxn
!)1(
)()( )1(
=
n
fxK xn?
=
-?=
n
i
i
x
n
n xxn
fxR
0
)1(
)(!)1( )()(?
上一页 下一页8
注,?通常不能确定?x,而是估计,?x?(a,b)
将 作为误差估计上限。
1)1( )( nn Mxf
=
-? n
i
in xxn
M
0
1 ||)!1(
当 f(x) 为任一个次数? n 的 多项式 时,,
可知,即插值多项式对于次数?n 的 多项式是 精确 的。
0)()1( xf n
0)(?xR n
Quiz:给定 xi = i +1,i = 0,1,2,3,4,5,下面哪个是 l2(x)的图像?
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
A B C
上一页 下一页9
例,已知
2 33s i n,214s i n,216s i n ===
分别利用 sin x 的 1次,2次 Lagrange 插值计算 sin 50?
并估计误差。
解:
0x 1x 2x
185500?=
n = 1 分别利用 x0,x1 以及 x1,x2 计算
4,6 10
== xx?利用
216/4/ 6/214/6/ 4/)(1?----= xxxL
这里 )
3,6(,s in)(,s in)( )2(-== xxxfxxf
而
)4)(6(!2 )()(,2 3s i n21 )2(1 --= xxfxR xx
0 0 7 6 2.0)185(0 1 3 1 9.0 1 --?R?sin 50?= 0.7660444…
)185(50sin 10L 0.77614
外推 /* extrapolation */ 的实际误差?-0.01010
3,4 21 == xx
利用 sin 50 0.76008,00660.0
185~00538.0 1R
内插 /* interpolation */ 的实际误差?0.00596
内插通常优于外推。选择要计算的 x 所在的区间的端点,插值效果较好。
上一页 下一页10
n = 2
2
3
))((
))((
2
1
))((
))((
21))((
))(()(
4363
46
3464
36
3646
342?-- ---- ---- --=
xxxxxxxL
)185(50sin 20L 0.76543
2 3c o s21;)3)(4)(6(!3
c o s)(
2---
-=
xx xxxxR
00 0 77.018500 0 44.0 2R?sin 50?= 0.7660444…
2次插值的实际误差?0.00061
高次插值通常优于低次插值但绝对不是次数越高就越好,……
上一页 下一页11
§ 2 牛顿插值 /* Newton’s Interpolation */
Lagrange 插值虽然易算,但若要增加一个节点时,
全部基函数 li(x) 都需重新算过。
将 Ln(x) 改写成,..))(()( 102010?--?-? xxxxaxxaa
)),,,(( 10 ---? nn xxxxa 的形式,希望每加一个节点时,
只附加一项 上去即可。
差商 (亦称均差 ) /* divided difference */
),()()(],[ ji
ji
ji
ji xxjixx
xfxfxxf
-
-=
1阶差商 /* the 1st
divided difference of f
w.r.t,xi and xj */
)(],[],[],,[ kixx xxfxxfxxxf
ki
kjji
kji?-
-= 2阶差商上一页 下一页12
1
11010
10
1110
10
],,...,[],,...,[
],,...,[],...,,[],...,[
--
-
-=
-
-=
kk
kkkk
k
kkk
k
xx
xxxfxxxf
xx
xxxfxxxfxxf
(k+1)阶差商:
=
= k
i ik
i
k x
xfxxf
0 1
0 )(
)(],.,,,[
事实上其中
,)()(
0
1?
=
-=
k
i
ik xxx
=
-=?
k
ij
j
jiik xxx
0
1 )()(?
差商的值与 xi 的顺序无关!
上一页 下一页13
牛顿插值 /* Newton’s Interpolation */
],[)()()( 000 xxfxxxfxf -?=
],,[)(],[],[ 101100 xxxfxxxxfxxf -?=
],...,,[)(],...,[],...,,[ 0010 nnnn xxxfxxxxfxxxf -?=-
)),,,((...))(()()( 10102010 -----?-?= nnn xxxxaxxxxaxxaaxN
1
2
… … … …
n-1
1 + (x - x0)? 2 + … … + ( x - x0)…( x - xn-1)? n-1
...))(](,,[)](,[)()( 102100100?--?-?= xxxxxxxfxxxxfxfxf
)),,,(](,...,[ 100 ---? nn xxxxxxf
))(),,,(](,...,,[ 100 nnn xxxxxxxxxf ---? -
Nn(x) R
n(x)
ai = f [ x0,…,xi ]
上一页 下一页14
注,?由 唯一性可知 Nn(x)? Ln(x),只是算法不同,故其余项也相同,即
)(!)1( )()(],...,,[ 1
)1(
10 xn
fxxxxf
k
x
n
kn?
=?
),(,! )(],...,[ m a xm i n)(0 xxkfxxf kk?=
实际计算过程为
f (x0)
f (x1)
f (x2)
…
f (xn-1)
f (xn)
f [x0,x1]
f [x1,x2]
… …
… …
f [xn-1,xn]
f [x0,x1,x2]
… …
… …
f [xn-2,xn-1,xn] f [x0,…,xn]
f (xn+1) f [xn,xn+1] f [xn-1,xn,xn+1] f [x1,…,xn+1] f [x0,…,xn+1]
上一页 下一页15
等距节点公式 /* Formulae with Equal Spacing */
向前差分
/* forward
difference */
iii fff -=1
i
k
i
k
i
k
i
k ffff 1
1
11 )( -
--?-?==?
向后差分
/* backward
difference */ 111 ---?-?=? ikikik fff
i-1ii fff -=?
中心差分
/* centered
difference */
2121
11 --?- -= ikikik fff
其中 )(
221 hii xff?=?
当节点 等距 分布时,),...,0(0 nihixx i =?=
上一页 下一页16
差分的重要性质:
线性:例如
gbfaxgbxfa= ))()((
若 f (x)是 m 次多项式,则 是 次多项式,而
)0()( mkxfk km-
)(0)( mkxfk?=?
差分值可由函数值算出:
=
-?-=?
n
j
jkn
j
k
n f
j
nf
0
)1(?
= -?
--=? n
j njk
jn
k
n f
j
nf
0
)1(
!
)1),,,(1(
j
jnnn
j
n?--=
其中 /* binomial coefficients */
函数值可由差分值算出,kj
n
j
kn fj
nf
=
=
0
k
k
k hk
fxxf
!],...,[
0
0
=
k
n
k
knnn hk
fxxxf
!],...,,[ 1
=
-- k
k
k
h
ff 0)( )(?=?
由 Rn 表达式上一页 下一页17
牛顿公式 )),,,(](,...,[...)](,[)()( 1000100 ----?= nnn xxxxxxfxxxxfxfxN
牛顿前差公式 /* Newton’s forward-difference formula */
牛顿后差公式 /* Newton’s backward-difference formula */
将节点顺序倒臵:
)),,,(](,...,[...)](,[)()( 101 xxxxxxfxxxxfxfxN nnnnnnn ---?= -
设 htxx?=
0
,则 )()()( 0
0
0 xfk
thtxNxN kn
k
nn?=?=?
=
),(,)),,,(1()!1( )()( 01)1( nnnn xxhntttnfxR?--?=
设 htxx
n?=
,则 )()1()()(
0
n
k
n
k
k
nnn xfk
thtxNxN?--=?=?
=
注,一般当 x 靠近 x0 时用前插,靠近 xn 时用后插,故两种公式亦称为 表初公式 和 表末公式 。
上一页 下一页18
§ 3 埃尔米特插值 /* Hermite Interpolation */
不仅要求函数值重合,而且要求若干阶 导数 也重合。
即:要求插值函数? (x) 满足?(xi) = f (xi),?’ (xi) = f ’ (xi),
…,?(mi) (xi) = f (mi) (xi).
注,?N 个条件可以确定 阶多项式。N - 1
要求在 1个节点 x0 处直到 m0阶导数都重合的插值多项式即为 Taylor多项式
0
0 )(
!
)(...))(()()(
00 0
)(
000
mm xxm xfxxxfxfx --=?
其余项为
)1(
0
0
)1(
0
0 )(
)!1(
)()()()( -
=-=
mm xx
m
fxxfxR
一般只考虑 f 与 f ’的值。
上一页 下一页19
例,设 x0? x1? x2,已知 f(x0),f(x1),f(x2) 和 f ’(x1),求多项式 P(x)
满足 P(xi) = f (xi),i = 0,1,2,且 P’(x1) = f ’(x1),并估计误差。
模仿 Lagrange 多项式的思想,设解,首先,P的阶数 = 3
= 2 13 )()()()()(
=0i ii
xhx1f ’xhxfxP?
h0(x) 有根 x1,x2,且 h0’(x1) = 0? x1 是重根。 )()()( 2
2
100 xxxxCxh --=
又,h0(x0) = 1? C0
)()(
)()()(
20
2
10
2
2
10
xxxx
xxxxxh
--
--=
h2(x)
h1(x) 有根 x0,x2? ))()(()( 201 xxxxBAxxh --?=
由余下条件 h1(x1) = 1 和 h1’(x1) = 0 可解。
与 h0(x) 完全类似。
(x)?h1 有根 x0,x1,x2h
1 ))()(()( 2101 xxxxxxCx ---=?
h1又,’(x1) = 1? C1 可解。
其中 hi(xj) =?ij,hi’(x1) = 0,(xi) = 0,’(x1) = 1?h1?h1
),())()(()()()( 221033 xxxxxxxKxPxfxR ---=-= !4 )()( )4(fxK?=
与 Lagrange 分析完全类似上一页 下一页20
一般地,已知 x0,…,xn 处有 y0,…,yn 和 y0’,…,yn’,求 H2n+1(x)
满足 H2n+1(xi) = yi,H’2n+1(xi) = yi’。
解,设= n i )()()(
=0i
i xhxhyixH2n+1
n
=0i
yi’
其中 hi(xj) =?ij,hi’(xj) = 0,(xj) = 0,’(xj) =?ij?hi?hi
hi(x) 有根 x0,…,xi,…,xn且都是 2重根 )()()( 2 xlBxAxh iiii?=
-
-=
ij ji
j
i xx
xxxl
)(
)()(
由余下条件 hi(xi) = 1 和 hi’(xi) = 0 可解 Ai 和 Bi?
)()])((21[)( 2 xlxxxlxh iiiii -?-=
(x)?hi 有根 x0,…,xn,除了 xi外都是 2重根hi )()( ii li2(x)xxCx -=
h
i又,’(xi) = 1? Ci = 1
h
i )(x )( i li2(x)xx-=
设 ],[,..,210 baCfbxxxa nn?== 则 2
0
)22(
)()!22( )()( -?=?
=
n
i
i
x
n
n xxn
fxR?
这样的 Hermite 插值唯一上一页 下一页21
Quiz,给定 xi = i +1,i = 0,1,2,3,4,5,下面哪个是 h2(x)的图像?
x0
-
-1
0.5
1 2 3 4 5 6
y
x
y
0
-
-
-1
0.5
1 2 3 4 5 6
斜率 =1?
求 Hermite多项式的基本步骤:
写出相应于条件的 hi(x),hi(x) 的组合形式;?
对每一个 hi(x),hi(x) 找出尽可能多的条件给出的根;?
根据多项式的总阶数和根的个数写出表达式;
根据尚未利用的条件解出表达式中的待定系数;
最后完整写出 H(x)。
上一页 下一页22
§ 4 分段低次插值 /* piecewise polynomial approximation */
例,在 [-5,5]上考察 的 Ln(x)。取
21
1)(
xxf?= ),.,,,0(105 niinx i =?-=
-5 -4 -3 -2 -1 0 1 2 3 4 5-0.5
0
0.5
1
1.5
2
2.5
n 越大,
端点附近抖动越大,称为
Runge 现象
Ln(x)? f (x)?
分段 低次 插值上一页 下一页23
分段线性插值 /* piecewise linear interpolation */
在每个区间 上,用 1阶多项式 (直线 ) 逼近 f (x):],[ 1?ii xx
11111 )()( -
-?
-
-=?
iii iiii i yxx
xxy
xx
xxxPxf ],[f o r
1 ii xxx
记,易证:当 时,||m a x
1 ii xxh -=? 0?h )()(1 xfxP h?
一致失去了原函数的光滑性。
上一页 下一页24
分段 Hermite插值 /* Hermite piecewise polynomials */
给定
nnn yyyyxx,...,;,...,;,...,000
在 上利用两点的 y 及 y’ 构造 3次 Hermite函数],[ 1?ii xx
导数一般不易得到。
上一页 下一页25
§ 5 三次样条 /* Cubic Spline */
定义 设 。 三次样条函数,
且在每个 上为 三次多项式 /* cubic polynomial */。 若它同时还满足,则称为 f 的 三次样条插值函数 /* cubic spline interpolant */.
bxxxa n ==,..10 ],[)( 2 baCxS?
],[ 1?ii xx
),...,0(),()( nixfxS ii ==
注,三次样条与分段 Hermite 插值的根本区别在于 S(x)自身光滑,不需要知道 f 的导数值(除了在 2个端点可能需要);而 Hermite插值依赖于 f 在所有插值点的导数值。
f(x) H(x)
S(x)
上一页 下一页26
构造三次样条插值函数的 三弯矩法
/* method of bending moment */
在 上,记],[
1 jj xx -,1--= jjj xxh ],[for )()( 1
][
jj
j xxxxSxS
-?=
对每个 j,此为 3次多项式则 S[j]”(x) 为 次多项式,需 个点的值确定之。1 2
设 S[j]”(xj-1) = Mj-1,S[j]”(xj) = Mj
对应力学中的 梁弯矩,故名对于 x?[xj-1,xj ] 可 得到
S[j]”(x) =
j
j
j
j
j
j h
xxM
h
xxM 1
1
-
-
-?-
积分 2次,可得 S[j]’(x) 和 S[j](x),
j
j
j
j
j
j
j Ah
xxM
h
xxM?-?-- -
-- 2
)(
2
)( 21
1
2
1S
[j]’(x) =
jj
j
j
j
j
j
j BxAh
xxM
h
xxM-?- -
- 6
)(
6
)( 313
1S
[j](x) =
利用已知
S[j](xj-1) = yj-1
S[j](xj) = yj
可解上一页 下一页27
j
jj
j
jj
j h
MM
h
yyA
6
11 -- ---=
j
j
j
j
j
j
j
j
j
jjj h
xxhMy
h
xxhMyBxA 1221
1 )6()6(
--
-
--?--=?
下面解决 Mj,利用 S’在 xj 的 连续性
[xj-1,xj ],S[j]’(x) = jjjjj
j
j
j
j
j
j h
MMxxf
h
xxM
h
xxM
6],[2
)(
2
)( 1
1
2
1
2
1
-
-
-
-
--?-?--
1
1
1
1
2
1
1
2
1
6],[2
)(
2
)(
--?-?--
j
jj
jj
j
j
j
j
j
j h
MMxxf
h
xxM
h
xxM[x
j,xj+1],S[j+1]’(x) =
利用 S[j]’(xj) = S[j+1]’(xj),合并关于 Mj-1,Mj,Mj+1的同类项,并记,,,整理后得到:
1
1
jj
j
j hh
h
=l 1 jj -= lm ]),[],[(
6
11
1
jjjj
jj
j xxfxxfhhg -?
-?=
2 11 gMMM jjjjjj =- lm
j?1 n-1
即:有 个未知数,个方程。n-1n+1
=
--- 1
1
0
11
11
2
2
n
n
nn
g
g
M
M
lm
lm
还需 2个 边界条件 /* boundary conditions */
上一页 下一页28
第 1类边条件 /* clamped boundary */,S’(a) = y0’,S’(b) = yn’
[a,x1 ],S[1]’(x) = 10110
1
2
1
1
2
1
0 6],[2
)(
2
)( hMMxxf
h
axM
h
xxM --?-?--
010
1
10 )],[(
62 gy
0’xxfhMM =-=?
nnn
n
nn gxxfyn’hMM =-=? -- ]),[(
62
11
类似地利用 [ xn-1,b ]
上的 S[n]’(x)
第 2类边条件,S”(a) = y0” = M0,S”(b) = yn” = Mn
这时:
nnn ygyg==== 2,0;2,0 000 ml
特别地,M0 = Mn = 0 称为 自由边界 /* free boundary */,对应的样条函数称为 自然样条 /* Natural Spline */。
第 3类边条件 /* periodic boundary */,
当 f 为 周期函数 时,
yn = y0,S’(a+) = S’(b-)
M0 = Mn
=
--
nnnn
nn
g
g
M
M 11
11
22
11
2
2
2
2
ml
lm
lm
ml
上一页 下一页29
注,?另有 三转角法 得到样条函数,即设 S[j]’(xj) = mj,则易知 [xj-1,xj ]上的 S[j](x) 就是 Hermite函数 。再利用 S”
的连续性,可导出关于 mj 的方程组,加上边界条件即可解。
Cubic Spline 由 boundary conditions 唯一 确定。
收敛性,若,且,则],[ baCf C
hh iim i nm a x一致
S(x)?f(x) 0m a x?
ihas
即,提高精度只须 增加节点,而无须提高样条阶数。
稳定性,只要边条件保证 | m0 |,| l0 |,| mn |,| l n | < 2,
则方程组系数阵为 SDD阵,保证数值稳定。
上一页 下一页30
Sketch of the Algorithm,Cubic Spline
① 计算 mj,l j,gj ; ② 计算 Mj (追赶法等 ) ;
③ 找到 x 所在区间 ( 即找到相应的 j ) ;
④ 由该区间上的 S[j](x) 算出 f(x) 的近似值。
插值法小结
Lagrange,给出 y0 … yn,选基函数 li(x),其次数为节点数 –1。
Newton? Ln(x),只是形式不同;节点等距或渐增节点时方便处理。
Hermite,给出 yi 及 yi ’,选 hi(x) 及 hi(x) 。
Spline:分段低次,自身光滑,f 的导数只在边界给出。
上一页 下一页31
作业,p.56-59
1题 (1),2题,4题,
5题,7题,9题,
10题,13题,20题。
上一页 下一页32
当精确函数 y = f(x) 非常复杂或未知时,在一系列节点 x0 … xn 处测得函数值 y0 = f(x0),…
yn = f(xn),由此构造一个简单易算的近似函数 g(x)? f(x),满足条件 g(xi) = f(xi) (i = 0,…
n)。这里的 g(x)称为 f(x) 的 插值函数 。最常用的插值函数是 …?多项式
x0 x1 x2 x3 x4x
g(x)? f(x)
上一页 下一页33
§ 1 拉格朗日多项式 /* Lagrange Polynomial */
niyxP iin,...,0,)( ==
求 n 次多项式 使得 n
nn xaxaaxP=?10)(
条件,无重合节点,即 ji xx?ji?
n = 1 已知 x0,x1 ; y0,y1,求 xaaxP
101 )(?=
使得
111001 )(,)( yxPyxP ==
可见 P1(x) 是过 ( x0,y0 ) 和 ( x1,y1 ) 两点的直线。
)()( 0
01
01
01 xxxx
yyyxP -
-
-?=
10
1
xx
xx
-
-
01
0
xx
xx
-
-= y
0 + y1
l0(x) l1(x)
== 1
0
)(
i ii
yxl
称为 拉氏基函数 /* Lagrange Basis */,
满足条件 li(xj)=?ij /* Kronecker Delta */
上一页 下一页34
n? 1 希望找到 li(x),i = 0,…,n 使得 li(xj)=?ij ;然后令
=
=
n
i
iin yxlxP
0
)()(,则显然有 Pn(xi) = yi 。
li(x) 每个 li 有 n 个根 x0 … xi … xn
=
-=---= n
j
j? i
jiniii xxCxxxxxxCxl
0
0 )())...()...(()(
-== j? i
ji
iii xxCxl )( 11)(
=
-
-= n
j
ij ji
j
i xx
xxxl
0
)(
)()(
=
= n
i
iin yxlxL
0
)()(
Lagrange
Polynomial
与 有关,而与 无关节点 f
上一页 下一页35
定理 (唯一性 ) 满足 的 n 阶插值多项式是唯一存在的。
niyxP ii,...,0,)( ==
证明,( p 25 利用 Vandermonde 行列式 论证 )
反证:若不唯一,则除了 Ln(x) 外还有另一 n 阶多项式 Pn(x) 满足 Pn(xi) = yi 。
考察 则 Qn 的阶数,)()()( xLxPxQ
nnn -=
n
而 Qn 有 个不同的根n + 1 x0 … xn
注,若不将多项式次数限制为 n,则插值多项式 不唯一 。
例如 也是一个插值多项式,其中 可以是任意多项式。
=
-?= n
i
in xxxpxLxP
0
)()()()(
)(xp
上一页 下一页36
插值余项 /* Remainder */
设节点
)1(?nf 在 [a,b]内存在,考察截断误差
)()()( xLxfxR nn -=
],[ baCf n?bxxxa
n10
,且 f 满足条件,
Rolle’s Theorem,若 充分光滑,,则存在 使得 。
)(x? 0)()( 10 == xx
),( 10 xx 0)( =
推广,若 0)()()(
210 === xxx ),(),,( 211100 xxxx
使得 0)()(
10 =?= ),( 10
使得 0)( =
0)()( 0 === nxx
存在 ),( ba 使得 0)()( = n
Rn(x) 至少有 个根n+1?
=
-=
n
i
in xxxKxR
0
)()()(
任意固定 x? xi (i = 0,…,n),考察?
=
--=
n
i
ixtxKtRnt
0
)()()()(?
(x)有 n+2 个不同的根 x0… xn x ),(,0)()1( ba
xxn?=
!)1()()()1(?-? nxKR xnn?
注意这里是对 t 求导
=?-- !)1)(()()( )1()1( nxKLf xnnxn
!)1(
)()( )1(
=
n
fxK xn?
=
-?=
n
i
i
x
n
n xxn
fxR
0
)1(
)(!)1( )()(?
上一页 下一页37
注,?通常不能确定?x,而是估计,?x?(a,b)
将 作为误差估计上限。
1)1( )( nn Mxf
=
-? n
i
in xxn
M
0
1 ||)!1(
当 f(x) 为任一个次数? n 的 多项式 时,,
可知,即插值多项式对于次数?n 的 多项式是 精确 的。
0)()1( xf n
0)(?xR n
Quiz:给定 xi = i +1,i = 0,1,2,3,4,5,下面哪个是 l2(x)的图像?
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
y
0
-
-
-1
0.5
-0.5
1 2 3 4 5 6 x
A B C
上一页 下一页38
例,已知
2 33s i n,214s i n,216s i n ===
分别利用 sin x 的 1次,2次 Lagrange 插值计算 sin 50?
并估计误差。
解:
0x 1x 2x
185500?=
n = 1 分别利用 x0,x1 以及 x1,x2 计算
4,6 10
== xx?利用
216/4/ 6/214/6/ 4/)(1?----= xxxL
这里 )
3,6(,s in)(,s in)( )2(-== xxxfxxf
而
)4)(6(!2 )()(,2 3s i n21 )2(1 --= xxfxR xx
0 0 7 6 2.0)185(0 1 3 1 9.0 1 --?R?sin 50?= 0.7660444…
)185(50sin 10L 0.77614
外推 /* extrapolation */ 的实际误差?-0.01010
3,4 21 == xx
利用 sin 50 0.76008,00660.0
185~00538.0 1R
内插 /* interpolation */ 的实际误差?0.00596
内插通常优于外推。选择要计算的 x 所在的区间的端点,插值效果较好。
上一页 下一页39
n = 2
2
3
))((
))((
2
1
))((
))((
21))((
))(()(
4363
46
3464
36
3646
342?-- ---- ---- --=
xxxxxxxL
)185(50sin 20L 0.76543
2 3c o s21;)3)(4)(6(!3
c o s)(
2---
-=
xx xxxxR
00 0 77.018500 0 44.0 2R?sin 50?= 0.7660444…
2次插值的实际误差?0.00061
高次插值通常优于低次插值但绝对不是次数越高就越好,嘿嘿 ……
上一页 下一页40
§ 2 牛顿插值 /* Newton’s Interpolation */
Lagrange 插值虽然易算,但若要增加一个节点时,
全部基函数 li(x) 都需重新算过。
将 Ln(x) 改写成,..))(()( 102010?--?-? xxxxaxxaa
)),,,(( 10 ---? nn xxxxa 的形式,希望每加一个节点时,
只附加一项 上去即可。
差商 (亦称均差 ) /* divided difference */
),()()(],[ ji
ji
ji
ji xxjixx
xfxfxxf
-
-=
1阶差商 /* the 1st
divided difference of f
w.r.t,xi and xj */
)(],[],[],,[ kixx xxfxxfxxxf
ki
kjji
kji?-
-= 2阶差商上一页 下一页41
1
11010
10
1110
10
],,...,[],,...,[
],,...,[],...,,[],...,[
--
-
-=
-
-=
kk
kkkk
k
kkk
k
xx
xxxfxxxf
xx
xxxfxxxfxxf
(k+1)阶差商:
=
= k
i ik
i
k x
xfxxf
0 1
0 )(
)(],.,,,[
事实上其中
,)()(
0
1?
=
-=
k
i
ik xxx
=
-=?
k
ij
j
jiik xxx
0
1 )()(?
差商的值与 xi 的顺序无关!
上一页 下一页42
牛顿插值 /* Newton’s Interpolation */
],[)()()( 000 xxfxxxfxf -?=
],,[)(],[],[ 101100 xxxfxxxxfxxf -?=
],...,,[)(],...,[],...,,[ 0010 nnnn xxxfxxxxfxxxf -?=-
)),,,((...))(()()( 10102010 -----?-?= nnn xxxxaxxxxaxxaaxN
1
2
… … … …
n-1
1 + (x - x0)? 2 + … … + ( x - x0)…( x - xn-1)? n-1
...))(](,,[)](,[)()( 102100100?--?-?= xxxxxxxfxxxxfxfxf
)),,,(](,...,[ 100 ---? nn xxxxxxf
))(),,,(](,...,,[ 100 nnn xxxxxxxxxf ---? -
Nn(x) R
n(x)
ai = f [ x0,…,xi ]
上一页 下一页43
注,?由 唯一性可知 Nn(x)? Ln(x),只是算法不同,故其余项也相同,即
)(!)1( )()(],...,,[ 1
)1(
10 xn
fxxxxf
k
x
n
kn?
=?
),(,! )(],...,[ m a xm i n)(0 xxkfxxf kk?=
实际计算过程为
f (x0)
f (x1)
f (x2)
…
f (xn-1)
f (xn)
f [x0,x1]
f [x1,x2]
… …
… …
f [xn-1,xn]
f [x0,x1,x2]
… …
… …
f [xn-2,xn-1,xn] f [x0,…,xn]
f (xn+1) f [xn,xn+1] f [xn-1,xn,xn+1] f [x1,…,xn+1] f [x0,…,xn+1]
上一页 下一页44
等距节点公式 /* Formulae with Equal Spacing */
向前差分
/* forward
difference */
iii fff -=1
i
k
i
k
i
k
i
k ffff 1
1
11 )( -
--?-?==?
向后差分
/* backward
difference */ 111 ---?-?=? ikikik fff
i-1ii fff -=?
中心差分
/* centered
difference */
2121
11 --?- -= ikikik fff
其中 )(
221 hii xff?=?
当节点 等距 分布时,),...,0(0 nihixx i =?=
上一页 下一页45
差分的重要性质:
线性:例如
gbfaxgbxfa= ))()((
若 f (x)是 m 次多项式,则 是 次多项式,而
)0()( mkxfk km-
)(0)( mkxfk?=?
差分值可由函数值算出:
=
-?-=?
n
j
jkn
j
k
n f
j
nf
0
)1(?
= -?
--=? n
j njk
jn
k
n f
j
nf
0
)1(
!
)1),,,(1(
j
jnnn
j
n?--=
其中 /* binomial coefficients */
函数值可由差分值算出,kj
n
j
kn fj
nf
=
=
0
k
k
k hk
fxxf
!],...,[
0
0
=
k
n
k
knnn hk
fxxxf
!],...,,[ 1
=
-- k
k
k
h
ff 0)( )(?=?
由 Rn 表达式上一页 下一页46
牛顿公式 )),,,(](,...,[...)](,[)()( 1000100 ----?= nnn xxxxxxfxxxxfxfxN
牛顿前差公式 /* Newton’s forward-difference formula */
牛顿后差公式 /* Newton’s backward-difference formula */
将节点顺序倒臵:
)),,,(](,...,[...)](,[)()( 101 xxxxxxfxxxxfxfxN nnnnnnn ---?= -
设 htxx?=
0
,则 )()()( 0
0
0 xfk
thtxNxN kn
k
nn?=?=?
=
),(,)),,,(1()!1( )()( 01)1( nnnn xxhntttnfxR?--?=
设 htxx
n?=
,则 )()1()()(
0
n
k
n
k
k
nnn xfk
thtxNxN?--=?=?
=
注,一般当 x 靠近 x0 时用前插,靠近 xn 时用后插,故两种公式亦称为 表初公式 和 表末公式 。
上一页 下一页47
§ 3 埃尔米特插值 /* Hermite Interpolation */
不仅要求函数值重合,而且要求若干阶 导数 也重合。
即:要求插值函数? (x) 满足?(xi) = f (xi),?’ (xi) = f ’ (xi),
…,?(mi) (xi) = f (mi) (xi).
注,?N 个条件可以确定 阶多项式。N - 1
要求在 1个节点 x0 处直到 m0阶导数都重合的插值多项式即为 Taylor多项式
0
0 )(
!
)(...))(()()(
00 0
)(
000
mm xxm xfxxxfxfx --=?
其余项为
)1(
0
0
)1(
0
0 )(
)!1(
)()()()( -
=-=
mm xx
m
fxxfxR
一般只考虑 f 与 f ’的值。
上一页 下一页48
例,设 x0? x1? x2,已知 f(x0),f(x1),f(x2) 和 f ’(x1),求多项式 P(x)
满足 P(xi) = f (xi),i = 0,1,2,且 P’(x1) = f ’(x1),并估计误差。
模仿 Lagrange 多项式的思想,设解,首先,P的阶数 = 3
= 2 13 )()()()()(
=0i ii
xhx1f ’xhxfxP?
h0(x) 有根 x1,x2,且 h0’(x1) = 0? x1 是重根。 )()()( 2
2
100 xxxxCxh --=
又,h0(x0) = 1? C0
)()(
)()()(
20
2
10
2
2
10
xxxx
xxxxxh
--
--=
h2(x)
h1(x) 有根 x0,x2? ))()(()( 201 xxxxBAxxh --?=
由余下条件 h1(x1) = 1 和 h1’(x1) = 0 可解。
与 h0(x) 完全类似。
(x)?h1 有根 x0,x1,x2h
1 ))()(()( 2101 xxxxxxCx ---=?
h1又,’(x1) = 1? C1 可解。
其中 hi(xj) =?ij,hi’(x1) = 0,(xi) = 0,’(x1) = 1?h1?h1
),())()(()()()( 221033 xxxxxxxKxPxfxR ---=-= !4 )()( )4(fxK?=
与 Lagrange 分析完全类似上一页 下一页49
一般地,已知 x0,…,xn 处有 y0,…,yn 和 y0’,…,yn’,求 H2n+1(x)
满足 H2n+1(xi) = yi,H’2n+1(xi) = yi’。
解,设= n i )()()(
=0i
i xhxhyixH2n+1
n
=0i
yi’
其中 hi(xj) =?ij,hi’(xj) = 0,(xj) = 0,’(xj) =?ij?hi?hi
hi(x) 有根 x0,…,xi,…,xn且都是 2重根 )()()( 2 xlBxAxh iiii?=
-
-=
ij ji
j
i xx
xxxl
)(
)()(
由余下条件 hi(xi) = 1 和 hi’(xi) = 0 可解 Ai 和 Bi?
)()])((21[)( 2 xlxxxlxh iiiii -?-=
(x)?hi 有根 x0,…,xn,除了 xi外都是 2重根hi )()( ii li2(x)xxCx -=
h
i又,’(xi) = 1? Ci = 1
h
i )(x )( i li2(x)xx-=
设 ],[,..,210 baCfbxxxa nn?== 则 2
0
)22(
)()!22( )()( -?=?
=
n
i
i
x
n
n xxn
fxR?
这样的 Hermite 插值唯一上一页 下一页50
Quiz,给定 xi = i +1,i = 0,1,2,3,4,5,下面哪个是 h2(x)的图像?
x0
-
-1
0.5
1 2 3 4 5 6
y
x
y
0
-
-
-1
0.5
1 2 3 4 5 6
斜率 =1?
求 Hermite多项式的基本步骤:
写出相应于条件的 hi(x),hi(x) 的组合形式;?
对每一个 hi(x),hi(x) 找出尽可能多的条件给出的根;?
根据多项式的总阶数和根的个数写出表达式;
根据尚未利用的条件解出表达式中的待定系数;
最后完整写出 H(x)。
上一页 下一页51
§ 4 分段低次插值 /* piecewise polynomial approximation */
例,在 [-5,5]上考察 的 Ln(x)。取
21
1)(
xxf?= ),.,,,0(105 niinx i =?-=
-5 -4 -3 -2 -1 0 1 2 3 4 5-0.5
0
0.5
1
1.5
2
2.5
n 越大,
端点附近抖动越大,称为
Runge 现象
Ln(x)? f (x)?
分段 低次 插值上一页 下一页52
分段线性插值 /* piecewise linear interpolation */
在每个区间 上,用 1阶多项式 (直线 ) 逼近 f (x):],[ 1?ii xx
11111 )()( -
-?
-
-=?
iii iiii i yxx
xxy
xx
xxxPxf ],[f o r
1 ii xxx
记,易证:当 时,||m a x
1 ii xxh -=? 0?h )()(1 xfxP h?
一致失去了原函数的光滑性。
上一页 下一页53
分段 Hermite插值 /* Hermite piecewise polynomials */
给定
nnn yyyyxx,...,;,...,;,...,000
在 上利用两点的 y 及 y’ 构造 3次 Hermite函数],[ 1?ii xx
导数一般不易得到。
上一页 下一页54
§ 5 三次样条 /* Cubic Spline */
定义 设 。 三次样条函数,
且在每个 上为 三次多项式 /* cubic polynomial */。 若它同时还满足,则称为 f 的 三次样条插值函数 /* cubic spline interpolant */.
bxxxa n ==,..10 ],[)( 2 baCxS?
],[ 1?ii xx
),...,0(),()( nixfxS ii ==
注,三次样条与分段 Hermite 插值的根本区别在于 S(x)自身光滑,不需要知道 f 的导数值(除了在 2个端点可能需要);而 Hermite插值依赖于 f 在所有插值点的导数值。
f(x) H(x)
S(x)
上一页 下一页55
构造三次样条插值函数的 三弯矩法
/* method of bending moment */
在 上,记],[
1 jj xx -,1--= jjj xxh ],[for )()( 1
][
jj
j xxxxSxS
-?=
对每个 j,此为 3次多项式则 S[j]”(x) 为 次多项式,需 个点的值确定之。1 2
设 S[j]”(xj-1) = Mj-1,S[j]”(xj) = Mj
对应力学中的 梁弯矩,故名对于 x?[xj-1,xj ] 可 得到
S[j]”(x) =
j
j
j
j
j
j h
xxM
h
xxM 1
1
-
-
-?-
积分 2次,可得 S[j]’(x) 和 S[j](x),
j
j
j
j
j
j
j Ah
xxM
h
xxM?-?-- -
-- 2
)(
2
)( 21
1
2
1S
[j]’(x) =
jj
j
j
j
j
j
j BxAh
xxM
h
xxM-?- -
- 6
)(
6
)( 313
1S
[j](x) =
利用已知
S[j](xj-1) = yj-1
S[j](xj) = yj
可解上一页 下一页56
j
jj
j
jj
j h
MM
h
yyA
6
11 -- ---=
j
j
j
j
j
j
j
j
j
jjj h
xxhMy
h
xxhMyBxA 1221
1 )6()6(
--
-
--?--=?
下面解决 Mj,利用 S’在 xj 的 连续性
[xj-1,xj ],S[j]’(x) = jjjjj
j
j
j
j
j
j h
MMxxf
h
xxM
h
xxM
6],[2
)(
2
)( 1
1
2
1
2
1
-
-
-
-
--?-?--
1
1
1
1
2
1
1
2
1
6],[2
)(
2
)(
--?-?--
j
jj
jj
j
j
j
j
j
j h
MMxxf
h
xxM
h
xxM[x
j,xj+1],S[j+1]’(x) =
利用 S[j]’(xj) = S[j+1]’(xj),合并关于 Mj-1,Mj,Mj+1的同类项,并记,,,整理后得到:
1
1
jj
j
j hh
h
=l 1 jj -= lm ]),[],[(
6
11
1
jjjj
jj
j xxfxxfhhg -?
-?=
2 11 gMMM jjjjjj =- lm
j?1 n-1
即:有 个未知数,个方程。n-1n+1
=
--- 1
1
0
11
11
2
2
n
n
nn
g
g
M
M
lm
lm
还需 2个 边界条件 /* boundary conditions */
上一页 下一页57
第 1类边条件 /* clamped boundary */,S’(a) = y0’,S’(b) = yn’
[a,x1 ],S[1]’(x) = 10110
1
2
1
1
2
1
0 6],[2
)(
2
)( hMMxxf
h
axM
h
xxM --?-?--
010
1
10 )],[(
62 gy
0’xxfhMM =-=?
nnn
n
nn gxxfyn’hMM =-=? -- ]),[(
62
11
类似地利用 [ xn-1,b ]
上的 S[n]’(x)
第 2类边条件,S”(a) = y0” = M0,S”(b) = yn” = Mn
这时:
nnn ygyg==== 2,0;2,0 000 ml
特别地,M0 = Mn = 0 称为 自由边界 /* free boundary */,对应的样条函数称为 自然样条 /* Natural Spline */。
第 3类边条件 /* periodic boundary */,
当 f 为 周期函数 时,
yn = y0,S’(a+) = S’(b-)
M0 = Mn
=
--
nnnn
nn
g
g
M
M 11
11
22
11
2
2
2
2
ml
lm
lm
ml
上一页 下一页58
注,?另有 三转角法 得到样条函数,即设 S[j]’(xj) = mj,则易知 [xj-1,xj ]上的 S[j](x) 就是 Hermite函数 。再利用 S”
的连续性,可导出关于 mj 的方程组,加上边界条件即可解。
Cubic Spline 由 boundary conditions 唯一 确定。
收敛性,若,且,则],[ baCf C
hh iim i nm a x一致
S(x)?f(x) 0m a x?
ihas
即,提高精度只须 增加节点,而无须提高样条阶数。
稳定性,只要边条件保证 | m0 |,| l0 |,| mn |,| l n | < 2,
则方程组系数阵为 SDD阵,保证数值稳定。
上一页 下一页59
Sketch of the Algorithm,Cubic Spline
① 计算 mj,l j,gj ; ② 计算 Mj (追赶法等 ) ;
③ 找到 x 所在区间 ( 即找到相应的 j ) ;
④ 由该区间上的 S[j](x) 算出 f(x) 的近似值。
插值法小结
Lagrange,给出 y0 … yn,选基函数 li(x),其次数为节点数 –1。
Newton? Ln(x),只是形式不同;节点等距或渐增节点时方便处理。
Hermite,给出 yi 及 yi ’,选 hi(x) 及 hi(x) 。
Spline:分段低次,自身光滑,f 的导数只在边界给出。
上一页 下一页60
§ 5 Cubic Spline
Lab 1,Cubic Spline
Construct the cubic spline interpolant S for the function f,defined at
points,satisfying some given boundary conditions,Partition
an given interval into m equal-length subintervals,and approximate the
function values at the endpoints of these subintervals.
Input
There are several sets of inputs,For each set:
The 1st line contains an integer 20? n? 0 which is the number of interpolating
points - 1,n = -1 signals the end of file.
The 2nd line contains n+1 distinct real numbers,
The 3rd line contains n+1 real numbers,
The 4th line contains an integer Type and three real numbers s0,sn,and Fmax.
If Type = 1,the clamped boundary condition is given,and,
If Type = 2,the natural boundary condition is given,and,
No other values of Type will be given.
Fmax is the default value that S (x) will assume if x is out of the range [x0,xn].
The last line of a test case consists of two real numbers t0 and tm,and an integer m?
1,Here the interval is partitioned into m equal-lengh subintervals,
…,The numbers are separated by spaces and new lines.
nxxx,...,,10
)(,...),(),( 10 nxfxfxf
nxxx,..10
00 )( sxS =? nn sxS =? )(
00 )( sxS = nn sxS = )(
],[ 0 mtt ],[ 10 tt
],[ 1 mm tt -
上一页 下一页61
§ 5 Cubic Spline
Output (? represents a space)
For each test case,you are supposed to output the following information:
1,The set of coefficients of S(x) in the format:
where
for,
Each coefficient is to be printed as in C printf:
fprintf(outfile,"%12.8e? ",coefficient );
2,m+1 approximating values of f at in the format:
fprintf(outfile,"f(%12.8e)? =? %12.8e\n",t,S );
The outputs of two test cases must be seperated by a blank line.
nnnn dcba
dcba
dcba
..
.
..
.
..
.
..
.
2222
1111
31211][ )()()()()( --- -?-?-?== jjjjjjjj xxdxxcxxbaxSxS
],[ 1 jj xxx -?
mtt,..,,0
Sample Input
2
0.0? 1.0? 2.0
0.0? 1.0? 2.0
1? 1.0? 1.0? 0.0
0.0? 3.0? 2
2
–0.5? –0.25? 0.0
–0.0247500? 0.3349375? 1.1010000
2? 0.0? 0.0? 0.0
–1.0? 0.0? 4
–1
上一页 下一页62
§ 5 Cubic Spline
Sample Output
0.00000000e+000? 1.00000000e+000? 0.00000000e+000? 0.00000000e+000?
1.00000000e+000? 1.00000000e+000? 0.00000000e+000? 0.00000000e+000?
f(0.00000000e+000)? =? 0.00000000e+000
f(1.50000000e+000)? =? 1.50000000e+000
f(3.00000000e+000)? =? 0.00000000e+000
–2.47500000e–002? 1.03237500e+000? 0.00000000e+000? 6.50200000e+000?
3.34937500e–001? 2.25150000e+000? 4.87650000e+000? –6.50200000e+000?
f(–1.00000000e+000)? =? 0.00000000e+000
f(–7.50000000e–001)? =? 0.00000000e+000
f(–5.00000000e–001)? =? –2.47500000e–002
f(–2.50000000e–001)? =? 3.34937500e–001
f(0.00000000e+000)? =? 1.10100000e+000