计算方法大作业——三次样条插值

(6)
要得到 S(x)的表达式，就要确定 Mi 的值，对于不同的边界条件，求解 Mi 的方程组 不同，第一类、第二类和第三类边界条件的三弯矩方程组分别为公式(7)~(9)。
2 2
1
2
2 n 2
2
n 1
M 1 d1 1M 0 d2 M 2 n 2 M n 2 dn2 d M 2 M n 1 n 1 n 2 n
1 0.9 0.8 0.7 0.6
S(x)
0.5 0.4 0.3 0.2 0.1 0 -1 -0.5 0 0.5 1
x
图 10 样条插值函数和插值节点图像
图 11 命令窗口运行结果
11
2 三次样条插值
②构造牛顿插值多项式 N5(x)和 N10(x)： 经过计算可得：
N5 x
125 4 45 2 59 x x 104 26 104
使用求三次样条插值函数的 matlab 程序对该问题进行求解，求解过程及结果如图 8 和图 9 所示。图 8 的结果与课本 137 页给出的结果完全一致，说明该 matlab 程序对 三次样条插值函数求解具有准确性，图 9 为将上式与插值节点画到同一个坐标系中得 到的图像，可以看到曲线非常光滑。
图 8 命令窗口显示的运行结果
9
2 三次样条插值
70 60 50 40
S(x)
30 20 10 0 -3
-2
-1
0
1
2
3
4
x
图 9 求解得到的三次样条插值函数与插值节点的图像
最终得到的结果为：
x3 9 x 2 25 x 28 3 2 x 3x 19 x 26 S x 3 2 2 x 3x 19 x 26 3 2 5 x 60 x 208 x 163
8
计算方法上机报告
此完成所有数据的输入。继续按 Enter 键会出现提示“选择封闭方程组的边界条件: 第 一类边界条件输入 1，第二类边界条件输入 2，第三类边界条件输入 3。 ”根据已知情况 选择相应的边界条件，若为自然三次样条插值，则选 1，并将插值区间两端点的二阶导 数值设置为 0。输入完成之后按 Enter 开始求解，程序运行结束后命令窗口会显示要求 的三次样条插值函数，同时会出现该插值函数以及插值节点的图像，便于直接观察。 2.3 算例及计算结果 (1) 《数值分析》课本第 137 页的例题 4.6.1，已知函数 y=f(x)的数值如下表，求它 的自然三次样条插值函数。 xi yi -3 7 -1 11 0 26 3 56 4 29
2 三次样条插值
2 三次样条插值
2.1 算法原理及程序框图 设在区间[a, b]上给定 n+1 个节点 xi（a ≤ x0 < x1 < … < xn ≤ b）,在节点 xi 处的函数 值为 yi = f(xi) (i = 0,1,…,n)。若函数 S(x)满足以下三个条件： (1) 在每个子区间[xi-1, xi] (i = 0,1,…,n)上，S(x)是三次多项式； (2) S(xi) = yi (i = 0,1,…,n)； (3) 在区间[a, b]上，S(x)的二阶导数 S”(x)连续， 则称 S(x)为函数 yi = f(x) 在区间[a, b]上的三次样条插值函数。 由定义可知 S(x)共有 4n 个待定参数，根据条件(3)可得如下 3n-3 个方程，
M 2 d 2 n 1 M n 1 d n 1 2 M n d n
(8)
1
2
1 M 1 d1
(9)
2 n 1
2
n
其中：hi xi xi 1 ，i
S x
x x i
6hi
3
M i 1
x xi 1
6hi
3
x x hi2 M i yi 1 M i 1 i 6 hi
x xi 1 h2 yi i M i , xi 1 x xi . 6 hi
(7)
6
计算方法上机报告
2 1
2 2 n
1 2
1
2
2
2 n 1
2 1
M 0 d0 M d 1 1 M 2 d2 d n 1 n 1 M n 1 2 M n d n
坐标系上进行比较。 ①用 matlab 程序求解三次样条插值函数 S10(x)（注意输入的差值节点数应为 11） ， 求解结果如图 10 和图 11 所示。最终得到的三次样条插值多项式为式(10)。
0.34171x3 1.02513 x 2 1.11327 x 0.46831 3 2 0.89326 x 2.34885 x 2.17225 x 0.75071 0.83639 x3 2.24648 x 2 2.11082 x 0.73842 3 2 13.40826 x 17.33272 x 8.14532 x 1.54302 3 2 54.46942 x 23.39388 x 1.00000 S ( x) 3 2 54.46942 x 23.39388 x 1.00000 13.40826 x3 17.33272 x 2 8.14532 x 1.54302 3 2 0.83639 x 2.24648 x 2.11082 x 0.73842 0.89326 x3 2.34885 x 2 2.17225 x 0.75071 3 2 0.34171x 1.02513 x 1.11327 x 0.46831
N10 x
22757 10 5444 8 20216 6 17147 4 3725 2 x x x x x 1 103 11 53 139 221
将牛顿插值多项式 N5(x)和 N10(x)及三次样条插值函数 S10(x)分别与 f(x)的曲线画在 同一个坐标系上进行比较，如图 12。可以看出三次样条函数与原函数符合的非常好， 对于低次的牛顿插值多项式，与原函数的大致趋势相同，而高次的牛顿插值多项式由 于龙格现象的出现，与原函数之间相差比较大。
7
2 三次样条插值
图 7 求三次样条插值函数的程序框图
2.2 程序使用说明 求三次样条插值函数的 matlab 程序见附录 3，该程序可以求解在给定区间上，连 续函数或者列表函数的三次样条插值函数。运行程序按照命令窗口的提示输入相关变 量直至得到结果。 运行该程序出现提示“选择差值节点的类型：差值对象为连续函数请输入 1，差值 对象为列表函数请输入 2。 ”根据已知条件输入 1 或者 2。若要求连续函数的三次样条 插值函数则输入 1；要求列表函数的插值函数则输入 2。当输入 1 时会出现提示“请输 入连续函数的表达式：f(x) =”此时需要输入连续函数的合法表达式，输入完成按 Enter 键， 接着还需要手动输入插值区间的下限 a 上限 b 以及节点个数 n； 当输入 2 时则要求 输入插值节点的个数 n、 插值节点构成的向量 x = (x1, x2,…, xn)以及插值节点函数值构成 的向量 y = (y1, y2,…, yn)，其中 x 和 y 输入的格式分别为[x1; x2;…; xn]和[y1; y2;…; yn]，至
1 x 0.8 0.8 x 0.6 0.6 x 0.4 0.4 x 0.2 0.2 x 0 0 x 0.2 0.2 x 0.4 0.4 x 0.6 0.6 x 0.8 0.8 x 1
(10)
10
计算方法上机报告
S ( xi ) S ( xi ), ( xi ) S ( xi ), S S ( x ) S ( x ), i i i 1, 2, , n 1
(5)
再加上条件(2)的 n+1 个方程， 共有 4n-2 个方程，因此还需要增加两个条件才可以确定 出 4n 个待定参数。所增加的条件称为边界条件或端点条件。常用的三种边界条件为： ①已知 f (a), f (b) ；②已知 f (a) f (b) ；③已知 f(x)是以 T=b-a 为周期的周期函数。 设 S(x)在节点 xi 处的二阶导数值为 Mi (i = 0,1,…,n)， Mi 为待定参数。 根据公式(5)可 以得到 S(x)在区间[xi-1, xi]上的表达式为
hn n 1 n ， ， hn h1
dn
6 y1 yn yn yn1 。 hn h1 h1 hn
三次样条插值函数的关键是要解式(7)、(8)或(9)所示的方程组，将得到的 Mi 代入
公式(6)中可得到相应的插值函数。其计算流程如图 7 所示。
(2) 给定函数 f ( x)
3 x 1 1 x 0 0 x3 3 x 4
1 (1 x 1) 。取等距节点，构造牛顿插值多项式 N5(x) 1 25x 2 和 N10(x)及三次样条插值函数 S10(x)。分别将三种插值多项式与 f(x)的曲线画在同一个
图 12 三个不同插值函数的比较
ห้องสมุดไป่ตู้12
hi ，i 1 i ，di 6 f [ xi 1 , xi , xi 1 ] ，i 1, 2, hi hi 1
, n 1 。
dn 6 f [ xn1 , xn , xn ] ； 公式(8)中的 d0 6 f [ x0 , x0 , x1 ] ， 公式(9)中的 n