§ 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)
§ 5 Cubic Spline
? 构造三次样条插值函数的 三弯矩法
/* method of bending moment */
在 上,记 ],[
1 jj xx ?,1??? jjj xxh ],[ for ) ( ) ( 1
] [
j j
j x x x x S x S
? ? ?
对每个 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
x x M
h
x x M 1
1
?
?
? ? ?
积分 2次,可得 S[j]’(x) 和 S[j](x),
j
j
j
j
j
j
j A h
x x M
h
x x M ? ? ? ? ? ?
? ? 2
) (
2
) ( 2 1
1
2
1 S
[j]’(x) =
j j
j
j
j
j
j
j B x A h
x x M
h
x x M ? ? ? ? ? ?
? 6
) (
6
) ( 3 1 3
1 S
[j](x) =
利用已知
S[j](xj?1) = yj?1
S[j](xj) = yj
可解
§ 5 Cubic Spline
j
j j
j
j j
j h
M M
h
y y A
6
1 1 ? ? ? ? ? ?
j
j
j
j
j
j
j
j
j
j j j h
x x h M y
h
x x h M y B x A 1 2 2 1
1 ) 6 ( ) 6 (
? ?
?
? ? ? ? ? ? ?
下面解决 Mj, 利用 S’ 在 xj 的 连续性
[xj?1,xj ],S[j]’(x) = j j j j j
j
j
j
j
j
j h
M M x x f
h
x x M
h
x x M
6 ],[ 2
) (
2
) ( 1
1
2
1
2
1
?
?
?
?
? ? ? ? ? ? ?
1
1
1
1
2
1
1
2
1
6 ],[ 2
) (
2
) (
?
?
?
?
?
?
? ? ? ? ? ? ? ?
j
j j
j j
j
j
j
j
j
j h
M M x x f
h
x x M
h
x x M [x
j,xj+1],S[j+1]’(x) =
利用 S[j]’(xj) = S[j+1]’(xj),合并关于 Mj?1,Mj,Mj+1的同类项,并
记,,,整理
后得到,
1
1
j j
j
j h h
h
?
?
? ? l 1 j j ? ? l m ]),[ ],[ (
6
1 1
1
j j j j
j j
j x x f x x f h h g ? ?
?
? ? ?
2 1 1 g M M M j j j j j j ? ? ? ? ? l m
? j ? 1 n?1
即:有 个未知数,个方程。 n?1 n+1
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??? 1
1
0
11
11
2
2
n
n
nn
g
g
M
M
lm
lm
还需 2个 边界条件 /* boundary conditions */
§ 5 Cubic Spline
? 第 1类边条件 /* clamped boundary */,S’(a) = y0’,S’(b) = yn’
[a,x1 ],S[1]’(x) = 1 0 1 1 0
1
2
1
1
2
1
0 6 ],[ 2
) (
2
) ( h M M x x f
h
a x M
h
x x M ? ? ? ? ? ? ?
0 1 0
1
1 0 ) ],[ (
6 2 g y
0’ x x f h M M ? ? ? ?
n n n
n
n n g x x f yn’ h M M ? ? ? ? ? ? ]),[ (
6 2
1 1
类似地利用 [ 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
§ 5 Cubic Spline
注,?另有 三转角法 得到样条函数,即设 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阵, 保证数值稳定。
HW,p.131 #4 #5
§ 5 Cubic Spline
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 的导数只在边界给出。
?
§ 5 Cubic Spline
Lab 11,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 ?
§ 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
§ 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
§ 6 快速傅立叶变换 /* Fast Fourier Transform */
? 问题的背景 /* background */
傅立叶变换 —— 函数展开为 三角级数
设 f(x) 周期为 2?,在 [0,2? ] 上展开为三角级数,
其中 Cj 为复系数,,则实际计算时要取
级数的前 N 项,并要求在区间的 N 个等分点上与 f(x)重合。
??
?0
)(
j
jxi
j eC
? ? ? ?xjixje jxi s i nco s)( ??
即:给定 [0,2? ] 上 N 个等分点 上的函
数值,令 满足插值条件 。
)1,...,1,0(2 ??? NkkNx k ?
)( kk xff ? ??
?
? 1
0
)()( N
j
jxij eCxS kk fxS ?)(
N 个未知数 N 个方程 )1,...,1,0(1
1
0
2 ??? ? ?
?
??????? Njef
NC
N
k
kNji
kj
?
Discrete
Four er
Transform
)1,...,1,0(
1
0
2
??? ?
?
?
?????? NkeCf N
j
kNji
jk
?
Inverse of
DFT
总之要进行形如
的计算,其中 已知,
??
?
? 1
0
N
k
jk
kj WxC
? ?kx ???????? NieW ?2
§ 6 Fast Fourier Transform
? Fast Fourier Transform
??
?
? 1
0
N
k
jk
kj WxC
???????? NieW ?2快速计算 ( j = 0,1,…,N?1),其中
直接计算需 复数乘法 次 N 2 ? 降到 N · logN
由于 W 的周期性 W qN+s = W s,W k j 实际上只有
这 个不同的值。若 N 为偶数,则 W k j只有 个
不同值。
1 0,.,? N W W
N N / 2
先 合并同类项,再做乘法。
I’m gonna need
some magic here …
§ 6 Fast Fourier Transform
例,N = 23 = 8,计算, j = 0,1,2,3,4,5,6,7 ?
?
? 7
0k
kjkj WxC
技巧,将 k,j 先记为 二进制数 /* binary numbers */
? ? ? ? ? ? ?
? ? ? ? ? ? ?
) ( 2 2 2
) ( 2 2 2
0 1 2
0
0
1
1
2
2
0 1 2
0
0
1
1
2
2
j j j j j j j
k k k k k k k
k j W ) 2 2 2 ( ) 2 2 2 ( 0 0 1 1 2 2 0 0 1 1 2 2 ? ? ? ? ? ? ? ? ? ? ? ? k k k j j j W
) ( ) 2 2 2 ( ) 2 2 2 ( 0 1 2 0 1 0 2 1 3 2 1 2 0 3 1 4 2 2 k k k j k k k j k k k j W ? ? ? ? ? ? ? ? ? ? ? ? ?
) ( ) 0 ( ) 0 0 ( 0 1 2 0 0 1 1 0 2 k k k j k k j k j W ? ? ?
? ? ?
? ? ? ??
?
?
?
??
?
?
?
?
?
?
?
?
?
?
1
0
)00(
1
0
)0(
1
0
)(
)()(
0
02
1
011
2
0120
012012
k
kj
k
kkj
k
kkkj
kkkjjj WWWxC
)( 0120 kkkA
)( 0011 jkkA
)( 0102 jjkA
)( 0123 jjjA
2?3次乘

全部计算需要 2?3 ?8次乘法
一般地:取 N = 2p,
每个 Cj 用 2p 次乘法,共用
2N log2N 次乘法。
利用,还可以进
一步化简到 N(p?1)/2 次乘法。
010 )1(2 jj pW ???
HW,
p.138 #1