《数值分析》课程教学资源（PPT讲稿）第六章 插值 nterpolationη（2/2）

§5三次样条 Cubic Spline 定义设a=x<x;<.<xn=b。三次样条函数S(x)eC1a,b, 且在每个x,x止为三次多项式/ cubic polynomial。若它同 时还满足S(x2)=f(x),(i=0,…,n)则称为∫的三次样条插值函 * / cubic spline interpolant*. 注:三次样条与分段 Hermite插值的根本区别在于(x)自 身光滑,不需要知道f的导数值(除了在2个端点可能需 要);而 hermite插值依赖于∫在所有插值点的导数值。 f(r) H() s(x)

§5 三次样条 /* Cubic Spline */ 定 义 设 。三次样条函数 , 且在每个 上为三次多项式 /* cubic polynomial */。若它同 时还满足 ，则称为 f 的三次样条插值函 数 /* cubic spline interpolant */. a  x0  x1  ...  xn  b ( ) [ , ] 2 S x  C a b [ , ] i i1 x x S(x ) f (x ), (i 0, ... ,n) i  i  注：三次样条与分段 Hermite 插值的根本区别在于S(x)自 身光滑，不需要知道 f 的导数值（除了在2个端点可能需 要）；而Hermite插值依赖于f 在所有插值点的导数值。 f(x) H(x) S(x)

§5 Cubic Spline >构造三次样条插值函数的三弯矩法 method of bending moment * 在|x,xl上,记b=x1<对应力学中的梁弯矩,故名 则Sux)为1次多项式,需2个卢a喱定之。 设S“(xx 对每个,此为次多项式 对于x∈x,xl可得到 SUP(x)=M+M 积分2次,可得SW(x)和Sw(x): SW(r)=-M(-x)2 x;1)2 +M 2h 2h+ 利用已知 Su(xi-1=y so(x)=M (x;-x) +M 可解 6 (x-x-1+Ax+B 6h

§5 Cubic Spline  构造三次样条插值函数的三弯矩法 /* method of bending moment */ 在 [ x j1 , x j ]上，记 , h j  x j  x j 1 ( ) ( ) for [ , ] 1 [ ] j j j S x  S x x  x  x 对每个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 S 1 [j]’(x) = j j j j j j j j A x B h x x M h x x M        6 ( ) 6 ( ) 3 1 3 S 1 [j](x) = 利用已知 S[j](xj1) = yj1 S[j](xj) = yj 可解

§5 Cubic spline A3iVi-1M; -M h;A x+B,=(; x3+( M1,、x 6 6 下面解决M:利用S在x的连续性 Ix-x so(r)=-M M +flr -1,x I Mi -Mi-1 2 M-M C小SH(x)=-= (x-x) +M +1 几 +1 6 利用S(x)=S+(x),合并关于M、M、M1的同类项,并 记4= b+n1,H=1-1,8=八x,x形1-几x,xD,整理 后得到:H1M1+2M1+1Mm1=8 M ≤j≤n-1 11 g1 即:有m+1个未知数,m-1个方程。 H-1 g 还需2个边界条件/ boundary conditions

§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 A x B y 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 f x x 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 f x x h x x M h x x [xj , xj+1]: S[j+1]’(x) = M 利用S[j]’(xj) = S[j+1]’(xj)，合并关于Mj1、 Mj、 Mj+1的同类项，并 记 , , , 整理 后得到： 1 1 j j j j h h h    l  j 1 j m   l ( [ , ] [ , ]) 6 1 1 1 j j j j j j j f x x f x x h h g       m jM j1  2M j  l jM j1  g j 1  j  n1 即：有n+1个未知数，n1个方程。                      1 1 0 1 1 1 1 2 2 n n n n g g M M m l m l 还需2个边界条件 /* boundary conditions */

§5 Cubic spline a第1类边条件/ clamped boundary+:s(a)=y0’,S(b)=yn M-M la, x1 SIT(x)=-Mo +M 2h1 类似地利用xn1,b 2M0+M1=x(f1x0,x1l-y)=g0 上的sm(x) M+2M (yn’-f|xn-1,xnD)=gn a第2类边条件:s"(a)=y”=M,S"(b)=yn”=M 这时:=0,S1=2y;An=0,gn=2yn 特别地,M=Mn=0称为自由边界/ free boundary+,对应的 样条函数称为自然样条/ Natural Spline a第3类边条件/ periodic boundary+: 当∫为周期函数时,「2 HTM1「g1 A222() yn=yo, S(a)=S(b) Mo=M u-I 2 a 2M

§5 Cubic Spline  第1类边条件 /* clamped boundary */： S’(a) = y0 ’ ， S’(b) = yn ’ [a , x1 ]: S[1]’(x) = 1 1 0 0 1 1 2 1 1 2 1 0 6 [ , ] 2 ( ) 2 ( ) h M M f x x h x a M h x x M        0 1 0 1 0 1 ( [ , ] ) 6 2 y g 0 f x x ’ h M  M    n n n n n n yn ’ f x x g h M  M      ( [ , ]) 6 1 2 1 类似地利用[ xn1 , b ] 上的 S[n]’(x)  第2类边条件： S”(a) = y0 ” = M0， S”(b) = yn ” = Mn 这时： n n n  0 , g  2y ;  0 , g  2y l0 0 0 m 特别地，M0 = Mn = 0 称为自由边界 /* free boundary */，对应的 样条函数称为自然样条 /* Natural Spline */。  第3类边条件 /* periodic boundary */ ： 当 f 为周期函数时， yn = y0 ， S’(a+) = S’(b)  M0 = Mn                n n n n n n g g M M1 1 1 1 2 2 1 1 2 2 2 2 l m m l m l l m

§5 Cubic spline 注:另有三转角法得到样条函数,即设Sm(x)=m,则 易知[x,x上的w(x)就是 Hermite函数。再利用s” 的连续性,可导出关于m的方程组,加上边界条件 即可解。 L Cubic Spline由 boundary conditions唯一确定。 矿收敛性:若∈Ca,且mb≤Cf(x) as max;→>0 即:提高精度只须增加节点,而无须提高样条阶数 稳定性:只要边条件保证|A,,n1,孔n|<2, 则方程组系数阵为SDD阵,保证数值稳定。 HW:p131#4#5

§5 Cubic Spline 注：另有三转角法得到样条函数，即设 S[j]’(xj) = mj，则 易知[xj1 , xj]上的S[j](x) 就是Hermite函数。再利用S” 的连续性，可导出关于mj的方程组，加上边界条件 即可解。 Cubic Spline 由boundary conditions 唯一确定。 收敛性：若 f C[a,b]，且 h  C   ，则 h i i min max 一致 S(x)  f(x) max  0 hi as 即:提高精度只须增加节点, 而无须提高样条阶数。 稳定性：只要边条件保证 | m0 |, | l0 |, | mn |, | l n | < 2， 则方程组系数阵为SDD阵，保证数值稳定。 HW: p.131 #4 #5

§5 Cubic spline Sketch of the Algorithm: Cubic spline ①计算两,λj,8 ②计算M(追赶法等) ③找到x所在区间(即找到相应的j); ④由该区间上的SU(x)算出f(x)的近似值。 插值法小结 ◆ Lagrange:给出v…,选基函数lx),其次数为 节点数=。 NEwton≡L(x),只是形式不同;节点等距或渐增节点 时方便处理。 Hermite:给出y及,选h)及A) ◆ Spline:分段低次,自身光滑,厂的导数只在边界给出

§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 的导数只在边界给出。 

85 Cubic Spline Lab ll Cubic Spline Construct the cubic spline interpolant S for the function f, defined at points o 1. Here the interval [to, t is partitioned into m equal-lengh subintervals, [to, t, Itm-1,tmI. The numbers are separated by spaces and new lines

§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 2 nd line contains n+1 distinct real numbers . The 3 rd 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. n x , x , ... , x 0 1 ( ), ( ), ... , ( ) 0 1 n f x f x f x n x  x  ...  x 0 1 0 0 S(x )  s n n S(x )  s 0 0 S(x )  s n n S(x )  s [ , ] 0 m t t [ , ] 0 1 t t [ , ] m 1 m t t 

§5 Cubic spline Output(a 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 b1 d b, c2 d, where S(x)=S(x)=a; +b; (x-x-1)+c, (x-xi-)+d(x-x -1) -1 b. c d Each coefficient is to be printed as in C printf: fprintf(outfile, %128e", coefficient ) 2. m+l approximating values of f at to, . ,tm in the format: fprintf(outfile, " f(%12 8e)0=0%12.8en",t, S); The outputs of two test cases must be seperated by a blank line. 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 0.5■-0.25■0.0 0.0247500■0.3349375■11010000 2■0.0■0.0■0.0 1.0■0.0■4

§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. n n n dn a b c a b c d a b c d . . . . . . . . . . . .2 2 2 2 1 1 1 1 3 1 2 1 1 [ ] ( ) ( ) ( ) ( ) ( )   j  j  j  j  j  j  j j S x S x a b x x c x x d x x [ , ] x  x j1 x j m t , ..., t 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.00000000■=■0.000000000 f(1500000000■=■1.500000 f(3.00000000■=■0.0000000.000 2.47500000e002■1.03237500e+000■0.00000000e+000■6.50200000e+000■ 3.34937500e001■225150000e+000■4.87650000e+000■-6.50200000e+000■ f(-1.00000000■=■0.0000000 f(-7.50000■=■0.00000000 f(-500000001)■=■-2.47500000002 f(-250000001)■=■3.34937500e001 f(0.000000=■1.10100000c000

§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 傅立叶变换—函数展开为三角级数 设/(x)周期为,在10?m1L江为三角级数∑Ce, 其中G为复 总之要进行形型xWM笔时要取 级数的前 的计算,已e (3)重合。 即:给定[0,∠ 上的函 数值f=∫(x) m疋值条件S(x)∫。 Discrete C=∑几c (24 k=0 0“yN器器 个未知数 f=∑C,e (2 (k=0,1,…,N-1) nverse of DET j=0

§6 快速傅立叶变换 /* Fast Fourier Transform */  问题的背景 /* background */ 傅立叶变换 —— 函数展开为三角级数 设 f(x) 周期为2，在 [0, 2 ] 上展开为三角级数 ， 其中Cj 为复系数， ，则实际计算时要取 级数的前 N 项，并要求在区间的N 个等分点上与f(x)重合。   0 ( ) j i jx j C e e j x i j x i jx cos sin ( )   即：给定[0, 2 ] 上N 个等分点 上的函 数值 ，令 满足插值条件 。 ( 0,1, ..., 1) 2  k k  N  N xk  ( ) k xk f  f     1 0 ( ) ( ) N j i jx j S x C e k k S( x )  f N 个未知数 N 个方程 ( 0, 1, ... , 1) 1 1 0 2              f e j N N C N k k N i j j k  Discrete Fourier Transform ( 0, 1, ... , 1) 1 0 2             f C e k N N j k N i j k j  Inverse of DFT 总之要进行形如 的计算，其中 已知，     1 0 N k k j C j xk W xk          N i W e 2