§ 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
][
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
可解
§ 5 Cubic Spline
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 */
§ 5 Cubic Spline
第 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
§ 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只有 个不同值。
10,.,?NWW
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 */
)(222
)(222
012
0
0
1
1
2
2
012
0
0
1
1
2
2
jjjjjjj
kkkkkkk
kjW )222()222( 001122001122 kkkjjjW
)()222()222( 012010213212031422 kkkjkkkjkkkjW
)()0()00( 012001102 kkkjkkjkjW
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
定义 设 。 三次样条函数,
且在每个 上为 三次多项式 /* 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
][
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
可解
§ 5 Cubic Spline
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 */
§ 5 Cubic Spline
第 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
§ 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只有 个不同值。
10,.,?NWW
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 */
)(222
)(222
012
0
0
1
1
2
2
012
0
0
1
1
2
2
jjjjjjj
kkkkkkk
kjW )222()222( 001122001122 kkkjjjW
)()222()222( 012010213212031422 kkkjkkkjkkkjW
)()0()00( 012001102 kkkjkkjkjW
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
