三次样条插值作业题
第三章 插值法 三次样条插值
问题
分段低次插值
在处理实际问题时，总是希望将所得到的数据点用得越多越好。
最简单的方法是用直线将函数值点直接连接。
分段低次插值
基本思想：用分段低次多项式来代替单个多项式。
具体作法：(1) 把整个插值区间分割成多个小区间；
(2) 在每个小区间上作低次插值多项式；
(3) 将所有插值多项式拼接整一个多项式。
优点：公式简单、运算量小、稳定性好、收敛性…
缺点：节点处的导数不连续，失去原函数的光滑性。
三次样条函数
样条函数
由一些按照某种光滑条件分段拼接起来的多项式组成的函数。
最常用的样条函数为三次样条函数，即由三次多项式组成，满足处处有二阶连续导数。
定义设节点a =x 0< x 1 < …< x n -1 < x n =b ，若函数
在每个小区间[x i , x i +1 ]上是三次多项式，则称其为三次样条函数。
如果同时满足s (x i ) = f (x i ) (i = 0, 1, 2, …, n )，则称s (x ) 为f (x ) 在[a , b ]上的三次样条函数。
],[)(2b a C x s ∈
利用线性插值公式，即可得的表达式：
求导得：
即：
：第一类边界条件（缺省边界条件）。
matlab三次样条插值例题解析
文章标题：深度解析Matlab三次样条插值1. 前言在数学和工程领域中，插值是一种常见的数值分析技术，它可以用来估计不连续数据点之间的值。
而三次样条插值作为一种常用的插值方法，在Matlab中有着广泛的应用。
本文将从简单到复杂，由浅入深地解析Matlab中的三次样条插值方法，以便读者更深入地理解这一技术。
2. 三次样条插值概述三次样条插值是一种利用分段三次多项式对数据点进行插值的方法。
在Matlab中，可以使用spline函数来进行三次样条插值。
该函数需要输入数据点的x和y坐标，然后可以根据需要进行插值操作。
3. 三次样条插值的基本原理在进行三次样条插值时，首先需要对数据点进行分段处理，然后在每个分段上构造出一个三次多项式函数。
这些多项式函数需要满足一定的插值条件，如在数据点处函数值相等、一阶导数相等等。
通过这些条件，可以得到一个关于数据点的插值函数。
4. Matlab中的三次样条插值实现在Matlab中，可以使用spline函数来进行三次样条插值。
通过传入数据点的x和y坐标，可以得到一个关于x的插值函数。
spline函数也支持在已知插值函数上进行插值点的求值，这为用户提供了极大的灵活性。
5. 三次样条插值的适用范围和局限性虽然三次样条插值在许多情况下都能够得到较好的插值效果，但也存在一些局限性。
在数据点分布不均匀或有较大噪音的情况下，三次样条插值可能会出现较大的误差。
在实际应用中，需要根据具体情况选择合适的插值方法。
6. 个人观点和总结通过对Matlab中三次样条插值的深度解析，我深刻地理解了这一插值方法的原理和实现方式。
在实际工程应用中，我会根据数据点的情况选择合适的插值方法，以确保得到准确且可靠的结果。
我也意识到插值方法的局限性，这为我在实际工作中的决策提供了重要的参考。
通过以上深度解析，相信读者已经对Matlab中的三次样条插值有了更加全面、深刻和灵活的理解。
在实际应用中，希望读者能够根据具体情况选择合适的插值方法，以提高工作效率和准确性。
计算方法大作业1 克服Runge现象
x3
x2
x
1
S1 ( x)
-0.34685
0.2086
0.073964
0.038462
S2 (x)
S (xi 0 ) S x(i 0 )
S
'
(xi
0) S
xi' (
0 )i
S
'
'
x(i
0)S
xi' ' (
0)
1 ,n2, . . . , 1
（1）
这里共有了 3n-3 个条件，再加上条件（2）中的 n+1 个插值条件，共有 4n-2 个条件，
因此还需要 2 个方程才能确定 S (x) .通常可在区间[a, b]的端点 a x0,b xn 上各加一个边
dn1
1
2
Mn
dn
（6）
2 1
2
2
2
1 M1 d1
M2
d2
n 1
2
n
1
M
n
1
dn1
n
n 2 M n dn
由式（1）内点拼接条件，可得
i M i1 2M i i M i1 d j i 1, 2,..., n 1
（3） （4）
其中
i
hi 1 hi1
, hi
i
hi hi 1
三次样条插值作业题
例1 设)(x f 为定义在[0,3]上的函数，有下列函数值表：且2.0)('0=x f ，1)('3-=x f ，试求区间[0,3]上满足上述条件的三次样条插值函数)(x s本算法求解出的三次样条插值函数将写成三弯矩方程的形式：)()6()()6()(6)(6)(211123131j j jj j j jj j j j jj j jj x x h h M y x x h h M y x x h M x x h M x s --+--+-+-=+++++其中，方程中的系数jj h M 6，jj h M 61+，jj j j h h M y )6(2-，jjj j h h M y )6(211++-将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。
以下为Matlab 代码：%============================= % 本段代码解决作业题的例1 %============================= clear all clc% 自变量x 与因变量y ，两个边界条件的取值 IndVar = [0, 1, 2, 3]; DepVar = [0, 0.5, 2, 1.5];LeftBoun = 0.2;RightBoun = -1;% 区间长度向量，其各元素为自变量各段的长度h = zeros(1, length(IndVar) - 1);for i = 1 : length(IndVar) - 1h(i) = IndVar(i + 1) - IndVar(i);end% 为向量μ赋值mu = zeros(1, length(h));for i = 1 : length(mu) - 1mu(i) = h(i) / (h(i) + h(i + 1));endmu(i + 1) = 1;% 为向量λ赋值lambda = zeros(1, length(h));lambda(1) = 1;for i = 2 : length(lambda)lambda(i) = h(i) / (h(i - 1) + h(i)); end% 为向量d赋值d = zeros(1, length(h) + 1);d(1) = 6 * ( (DepVar(2) - DepVar(1) ) / ( IndVar(2) - IndVar(1) ) - LeftBoun) / h(1); for i = 2 : length(h)a = ( DepVar(i) - DepVar(i - 1) ) / ( IndVar(i) - IndVar(i - 1) );b = ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) );c = (b - a) / ( IndVar(i + 1) - IndVar(i - 1) );d(i) = 6 * c;endd(i + 1) = 6 *( RightBoun - ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) ) ) / h(i);% 为矩阵A赋值% 将主对角线上的元素全部置为2A = zeros( length(d), length(d) );for i = 1 : length(d)A(i, i) = 2;end% 将向量λ的各元素赋给主对角线右侧第一条对角线for i = 1 : length(d) - 1A(i, i + 1) = lambda(i);end% 将向量d的各元素赋给主对角线左侧第一条对角线for i = 1 : length(d) - 1A(i + 1, i) = mu(i);end% 求解向量MM =A \ d';% 求解每一段曲线的函数表达式for i = 1 : length(h)Coefs_1 = M(i) / (6 * h(i));Part_1 = conv( Coefs_1, ...conv( [-1, IndVar(i + 1)], ...conv( [-1, IndVar(i + 1)], [-1, IndVar(i + 1)] ) ) ); S_1 = polyval (Part_1, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_2 = M(i + 1)/(6 * h(i));Part_2 = conv( Coefs_2, ...conv( [1, -IndVar(i)], ...conv( [1, -IndVar(i)], [1, -IndVar(i)] ) ) );S_2 = polyval (Part_2, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_3 = (DepVar(i) - M(i) * h(i)^2 / 6) / h(i);Part_3 = conv(Coefs_3, [-1, IndVar(i + 1)]);S_3 = polyval (Part_3, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_4 = (DepVar(i + 1) - M(i + 1) * h(i)^2 / 6) / h(i);Part_4 = conv(Coefs_4, [1, -IndVar(i)]);S_4 = polyval (Part_4, [IndVar(i) : 0.01 : IndVar(i + 1)]);S = S_1 + S_2 + S_3 + S_4;plot ([IndVar(i) : 0.01 : IndVar(i + 1)], S, 'LineWidth', 1.25)% 在样条插值曲线的相应位置标注该段曲线的函数表达式text(i - 1, polyval(Part_1, 3), ...['\itS', num2str(i), '(x)=', num2str(Coefs_1), '(', num2str( IndVar(i + 1) ), '-x)^{3}+', ...num2str(Coefs_2), '(x-', num2str( IndVar(i) ), ')^{3}+', num2str(Coefs_3), ...'(', num2str( IndVar(i + 1) ), '-x)+', num2str(Coefs_4), '(x-',num2str( IndVar(i) ), ')'], ...'FontName', 'Times New Roman', 'FontSize', 14)hold onend% 过x=1和x=2两个横轴点作垂线 %line([1, 1], [2.5, -0.5], 'LineStyle', '--');line([2, 2], [2.5, -0.5], 'LineStyle', '--');% 为x轴和y轴添加标注xlabel( '\itx', 'FontName', 'Times New Roman', ...'FontSize', 14, 'FontWeight', 'bold');ylabel( '\its(x)', 'FontName', 'Times New Roman', ...'Rotation', 0, 'FontSize', 14, 'FontWeight', 'bold');最终，三次样条插值函数s(x)表达式为：[][][]⎪⎩⎪⎨⎧∈-+-+-+--∈-+-+---∈+-++--=.3,2,)2(44.1)3(62.2)2(06.0)3(62.0,2,1,)1(62.2)2(08.0)1(62.0)2(42.0,1,0,08.0)1(06.042.0)1(06.0)(333333x x x x x x x x x x x x x x x x s曲线的图像如图所示：例2 已知函数值表：试求在区间[1,5]上满足上述函数表所给出的插值条件的三次自然样条插值函数)(x s本算法求解出的三次样条插值函数将写成三弯矩方程的形式：)()6()()6()(6)(6)(211123131j j jj j j jj j j j jj j jj x x h h M y x x h h M y x x h M x x h M x s --+--+-+-=+++++其中，方程中的系数jj h M 6，jj h M 61+，jj j j h h M y )6(2-，jjj j h h M y )6(211++-将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。
三次样条插值的求解
三次样条插值的求解摘要：分段低次插值虽然解决了高次插值的振荡现象和数值不稳定现象，使得插值多项式具有一致收敛性，保证了插值函数整体的连续性，但在函数插值节点处不能很好地保证光滑性要求，这在某些要求光滑性的工程应用中是不能接受的。
如飞机的机翼一般要求使用流线形设计，以减少空气阻力，还有船体放样等的型值线，往往要求有二阶光滑度（即有二阶连续导数）。
因此，在分段插值的基础上，引进了一种新的插值方法，在保证原方法的收敛性和稳定性的同时，又使得函数具有较高的光滑性的样条插值。
关键字:三转角方程 三弯矩阵方程0. 引言1，三次样条函数定义1：若函数2()[,]S x a b C ∈,且在每个小区间上1,j j x x +⎡⎤⎦⎣上是三次多项式，其中01n a x x x b ⋯=<<<= 是给定节点，则称()s x 是节点01,,,n x x x ⋯上的三次样条函数。
若节点j x 上 给定函数值()j j y f x =(0,1,)j n ⋯= ,且()j j s x y = (0,1,)j n ⋯= （1.1）成立，则称 ()s x 为三次样条差值函数。
从定义知，要求出()s x ，在每个应小区间1[,]j j x x + 上确定4个待定系数，共有 n 个小区间，故应确定4n 个参数，根据()s x 在[,]a b 上二阶导数连续，在节点()1,2,3,,1j x j n ⋯=-处应满足连续性条件(0)(0),j j s x s x -=+ ''(0)(0),j j s x s x -=+''''(0)(0)j j s x s x -=+ （1.2） 共有 3n-3个条件，再加上()s x 满足插值条件（1.1），共有4n-2个条件，因此还需要2个条件才能确定()s x 。
通常可在区间[,]a b 端点0,n a x b x ==上各加一个条件（称边界条件），边界条件可根据实际的问题要求给定。
计算方法大作业——三次样条插值
计算方法上机报告
此完成所有数据的输入。继续按 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 个方程，
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
三次样条插值作业题
例1 设)(x f 为定义在[0,3]上的函数，有下列函数值表：且2.0)('0=x f ，1)('3-=x f ，试求区间[0,3]上满足上述条件的三次样条插值函数)(x s本算法求解出的三次样条插值函数将写成三弯矩方程的形式：)()6()()6()(6)(6)(211123131j j jj j j jj j j j jj j jj x x h h M y x x h h M y x x h M x x h M x s --+--+-+-=+++++其中，方程中的系数jj h M 6，jj h M 61+，jj j j h h M y )6(2-，jjj j h h M y )6(211++-将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。
以下为Matlab 代码：%=============================% 本段代码解决作业题的例1%============================= clear all clc% 自变量x 与因变量y ，两个边界条件的取值 IndVar = [0, 1, 2, 3]; DepVar = [0, 0.5, 2, 1.5];LeftBoun = 0.2; RightBoun = -1;% 区间长度向量，其各元素为自变量各段的长度 h = zeros(1, length(IndVar) - 1); for i = 1 : length(IndVar) - 1h(i) = IndVar(i + 1) - IndVar(i); end% 为向量μ赋值mu = zeros(1, length(h));for i = 1 : length(mu) - 1mu(i) = h(i) / (h(i) + h(i + 1));endmu(i + 1) = 1;% 为向量λ赋值lambda = zeros(1, length(h));lambda(1) = 1;for i = 2 : length(lambda)lambda(i) = h(i) / (h(i - 1) + h(i));end% 为向量d赋值d = zeros(1, length(h) + 1);d(1) = 6 * ( (DepV ar(2) - DepVar(1) ) / ( IndVar(2) - IndVar(1) ) - LeftBoun) / h(1);for i = 2 : length(h)a = ( DepVar(i) - DepVar(i - 1) ) / ( IndVar(i) - IndVar(i - 1) );b = ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) );c = (b - a) / ( IndVar(i + 1) - IndVar(i - 1) );d(i) = 6 * c;endd(i + 1) = 6 *( RightBoun - ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) ) ) / h(i);% 为矩阵A赋值% 将主对角线上的元素全部置为2A = zeros( length(d), length(d) );for i = 1 : length(d)A(i, i) = 2;end% 将向量λ的各元素赋给主对角线右侧第一条对角线for i = 1 : length(d) - 1A(i, i + 1) = lambda(i);end% 将向量d的各元素赋给主对角线左侧第一条对角线for i = 1 : length(d) - 1A(i + 1, i) = mu(i);end% 求解向量MM =A \ d';% 求解每一段曲线的函数表达式for i = 1 : length(h)Coefs_1 = M(i) / (6 * h(i));Part_1 = conv( Coefs_1, ...conv( [-1, IndVar(i + 1)], ...conv( [-1, IndVar(i + 1)], [-1, IndVar(i + 1)] ) ) );S_1 = polyval (Part_1, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_2 = M(i + 1)/(6 * h(i));Part_2 = conv( Coefs_2, ...conv( [1, -IndVar(i)], ...conv( [1, -IndVar(i)], [1, -IndVar(i)] ) ) );S_2 = polyval (Part_2, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_3 = (DepVar(i) - M(i) * h(i)^2 / 6) / h(i);Part_3 = conv(Coefs_3, [-1, IndVar(i + 1)]);S_3 = polyval (Part_3, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_4 = (DepVar(i + 1) - M(i + 1) * h(i)^2 / 6) / h(i);Part_4 = conv(Coefs_4, [1, -IndVar(i)]);S_4 = polyval (Part_4, [IndVar(i) : 0.01 : IndVar(i + 1)]);S = S_1 + S_2 + S_3 + S_4;plot ([IndVar(i) : 0.01 : IndVar(i + 1)], S, 'LineWidth', 1.25)% 在样条插值曲线的相应位置标注该段曲线的函数表达式text(i - 1, polyval(Part_1, 3), ...['\itS', num2str(i), '(x)=', num2str(Coefs_1), '(', num2str( IndVar(i + 1) ), '-x)^{3}+', ...num2str(Coefs_2), '(x-', num2str( IndVar(i) ), ')^{3}+', num2str(Coefs_3), ...'(', num2str( IndVar(i + 1) ), '-x)+', num2str(Coefs_4), '(x-', num2str( IndVar(i) ), ')'], ...'FontName', 'Times New Roman', 'FontSize', 14)hold onend% 过x=1和x=2两个横轴点作垂线%line([1, 1], [2.5, -0.5], 'LineStyle', '--');line([2, 2], [2.5, -0.5], 'LineStyle', '--');% 为x轴和y轴添加标注xlabel( '\itx', 'FontName', 'Times New Roman', ...'FontSize', 14, 'FontWeight', 'bold');ylabel( '\its(x)', 'FontName', 'Times New Roman', ...'Rotation', 0, 'FontSize', 14, 'FontWeight', 'bold');最终，三次样条插值函数s(x)表达式为：[][][]⎪⎩⎪⎨⎧∈-+-+-+--∈-+-+---∈+-++--=.3,2,)2(44.1)3(62.2)2(06.0)3(62.0,2,1,)1(62.2)2(08.0)1(62.0)2(42.0,1,0,08.0)1(06.042.0)1(06.0)(333333x x x x x x x x x x x x x x x x s曲线的图像如图所示：例2 已知函数值表：试求在区间[1,5]上满足上述函数表所给出的插值条件的三次自然样条插值函数)(x s本算法求解出的三次样条插值函数将写成三弯矩方程的形式：)()6()()6()(6)(6)(211123131j j jj j j jj j j j jj j jj x x h h M y x x h h M y x x h M x x h M x s --+--+-+-=+++++其中，方程中的系数jj h M 6，jj h M 61+，jj j j h h M y )6(2-，jjj j h h M y )6(211++-将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。
MATLAB大作业 给定一个时间序列,使用三次样条插值方法进行均匀内插
MATLAB作业给定一个时间序列，使用三次样条插值方法进行均匀内插（题目的相关说明：按题目要求编写一个MATLAB程序函数，并把自己编制程序所得的结果与MATLAB库函数分析结果进行对比。
）理论基础：时间序列的概念：时间序列是一种定量预测方法，又称简单外延法，时间序列分析是根据系统观测得到的时间序列数据，通过曲线拟合和参数估计来建立数学模型的理论与方法，时间序列分析可分为以下三种情况（1）把一个时间序列的值变动为N 个组成部分，通常可以分为四种 a、倾向变动，又称长期趋势变动 b、循环变动，又称周期变动 c、季节变动，即每年有规则的反复进行变动 d、不规则变动，即随机变动。
然后把这四个综合到一起得出预测的结果。
虽然分成这四部分，但这四部分之间的相互关系是怎么样的呢，目前一般采用相乘的关系，其实各个部分都是在其他部分作用的基础上进行作用的，所以采用相乘是有一定依据的，此种方法适合于短期预测和库存预测（2）把预测对象、预测目标和对预测的影响因素都看成为具有时序的，为时间的函数，而时间序列法就是研究预测对象自身变化过程及发展趋势，如果未来预测是线性的，其数学模型为YT+L=aT+bTL,YT+L为未来预测值，aT为截距，bT为斜率，L为由T到需要预测的单位时间数（如5年、10年等）（3）根据预测对象与影响因素之间的关系及影响程度来推算未来，与目标的相关因素很多，只能选择那些因果关系较强的为预测影响的因素，此时间序列法用于短期预测比较有效，若要用于长期预测，还需要结合其他方法才行。
三次样条插值的实际应用：在制造船体和汽车外形等工艺中传统的设计方法是，首先由设计人员按外形要求，给出外形曲线的一组离散点值，施工人员准备好有弹性的样条（一般用竹条或有弹性的钢条）和压铁，将压铁放在点的位置上，调整竹条的形状，使其自然光滑，这时竹条表示一条插值曲线，我们称为样条函数。
从数学上看，这一条近似于分段的三次多项式，在节点处具有一阶和二阶连续微商。
第五章(3)三次样条插值
6( xi xi 1 2 x ) ( yi 1 yi ) 3 hi 1
而
2 4 6 S ( xi 0) mi 1 mi 2 ( yi yi 1 ) hi hi hi 4 2 6 S ( xi 0) mi m i 1 2 ( yi 1 yi ) hi 1 hi 1 hi 1
n
当n 时，Ln ( x )只在 | x | 3.63 内收敛，而在该区间外 是发散的。
从图中可以看出，在 0 附近插值效果是好的，即余项较 小，另一种现象是插值多项式随节点增多而振动。这种插值 多项式当节点增加时反而不能更好地接近被插值函数的现象， 称为龙格现象。
上述现象告诉我们用高次插值多项式是不 妥当的，从数值计算上可解释为高次插值多项 式的计算会带来舍入误差的增大，从而引起计 算失真。因此，实践上作插值时一般只用一次、 二次最多用三次插值多项式。
式中x [ xi 1 , xi ] (i 1,2,, n)
第(2)步
为了确定mi，需要用到S ( x )的二阶导数在节点连续 的条件， S ( x )在[ xi 1 , xi ]和[ xi , xi 1 ]上的二阶导数分别为
Si( x ) 6 x 2 xi 1 4 xi 6 x 4 xi 1 2 xi mi 1 mi 2 2 hi hi ( x [ xi 1 , xi ])
若记hi xi xi 1，则上式可写为
( x x i ) 2 hi 2( x x i 1 ) ( x x i 1 ) 2 hi 2( x i x ) Si ( x) y i 1 yi 3 3 hi hi ( x x i ) 2 ( x x i 1 ) ( x x i 1 ) 2 ( x x i ) m i 1 mi 2 2 hi hi
_三次样条插值1
能出现Runge现象,采用分段插值虽然计算简单、也有
一致收敛性,但不能保证整条曲线在连接点处的光滑性 ,如分段线性插值,其图形是锯齿形的折线,虽然连续,但 处都是“尖点”,因而一阶导数都不存在,这在实用上, 往往不能满足某些工程技术的高精度要求。如在船体、
(5.39)
其中
6 g 0 h ( f x0 , x1 y 0 ) 1 g 6 ( y f x , x ) n n 1 n n hn
第二种边界条件:即已知插值区间两端的二阶导数值: S ( x0, y0, S ( xn ) yn ) 由于在区间端点处二阶导数
y n y n 1 6 ( yn ) hn hn
(5.38)
将式（5.36）和式（5.37）以及式（5.38）合在一起 即得确定 M 0 , M 1 ,, M n 的线性方程组
2 1 1 2
1
n 1
2 1
M 0 g0 M g 1 1 n 1 M n 1 g n 1 2 M n g n
其中,Ai,Bi为积分常数,可利用插值条件
S ( xi 1 ) f ( xi 1 ), S ( xi ) f ( xi ) 确定,即要求Ai,Bi满足
1 S ( xi 1 ) M i 1 hi2 Ai hi f ( xi 1 ), 6 1 S ( xi ) M i hi2 Bi hi f ( xi ) 6
在左端点xi-1上有
S i( xi 1 0) hi h y y h h y y M i 1 i (M i M i 1 ) i i 1 i M i 1 i M i i i 1 3 6 hi 2 6 hi
数值分析作业-三次样条插值
数值计算方法作业实验4.3 三次样条差值函数实验目的：掌握三次样条插值函数的三弯矩方法。
实验函数:dt ex f xt ⎰∞--=2221)(π实验容：（1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值；（3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线比较插值结果。
实验4.5 三次样条差值函数的收敛性实验目的：多项式插值不一定是收敛的，即插值的节点多，效果不一定好。
对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。
实验容：按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。
实验要求：（1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；（2） 三次样条插值函数的思想最早产生于工业部门。
作为工业应用的例子，考虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一算法描述：拉格朗日插值：其中是拉格朗日基函数，其表达式为：()∏≠=--=ni j j j i ji x x x x x l 0)()(牛顿插值：))...()(](,...,,[....))(0](,,[)0](,[)()(1102101210100----++--+-+=n n n x x x x x x x x x x f x x x x x x x f x x x x f x f x N其中⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧--=--=--=-)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[01102110x x x x x f x x x f x x x f x x x x f x x f x x x f x x x f x f x x f n n n n i k j i k j k j i ji j i j i三样条插值：所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a<X0<X1……<Xn<b)分成的每个小区间[x i-1,x i ]上是三次多项式，其在此区间上的表达式如下：],[),6()6(]6)([6)(6)()(111113131i i ii i i i i i i i i i i i i i i i i i x x x h yM h M h h y x M M h h y y h x x Mi h x x M x S -------∈-+-+---+-+-=式中Mi=)(i x S ''.因此，只要确定了Mi 的值，就确定了整个表达式，Mi 的计算方法如下：令⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=---+=+=+=+--++++++],,[6)(6111111111i i i i i i i i i i i i i i i i i i i ix x x f h y y h y y h h d h h h h h h λμ则Mi 满足如下n-1个方程：1,...2,1,211-==+++-n i d M M M i i i i i i λμ常用的边界条件有如下几类：（1） 给定区间两端点的斜率m 0,m n ,即n n n m y x S m y x S ='='='=')(,)(000 （2） 给定区间两端点的二阶导数M0，Mn,即n n n M y x S M y x S =''=''=''='')(,)(000 （3） 假设y=f(x)是以b-a 为周期的周期函数，则要求三次样条插值函数S （x ）也为周期函数，对S （x ）加上周期条件2,1,0),0()0()(0)(=-=+p x S x S n p p对于第一类边界条件有⎪⎪⎩⎪⎪⎨⎧--=+--=+--)(62)(6211001110n n n n n n ih y y mn h M M m hy y h M M对于第二类边界条件有⎩⎨⎧=+=+-n n n n d M M d M M 221100μλ其中n n n n nnn M u x x f m h d M m x x f h d )1(2]),[(6)1(2)],[(6100001010-+-=-+-=-μλλ那么解就可以为⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n n n n d d d d d M M M M M 1210121011...2...............2............................1..2.1......0..2μλμλμλ 对于第三类边界条件，)0()0(,,000-=+==n n n x S x S M M y y ，由此推得0010012d M M M n =-++μλ，其中]),1[],[(6,,101010110n n nn n n x x f x x f h h d h h h h h h --+=+=+=μλ，那么解就可以为： ⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-------1221012101221100...2.............2..............................2..,,.......,..22n n n n n n n d d d d d M M M M M n μλλμλμμλ 程序代码： 1拉格朗日插值函数Lang.mfunction f=lang(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标 %xi 插值点处的横坐标%f 求得的拉格朗日插值多项式的值 n=length(X); f=0; for i=1:n l=1; for j=1:i-1l=l.*(xi-X(j))/(X(i)-X(j)); end ; for j=i+1:nl=l.*(xi-X(j))/(X(i)-X(j)); end ;%拉格朗日基函数 f=f+l*Y(i); endfprintf('%d\n',f) return2 牛顿插值函数newton.mfunction f=newton(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标 %xi 插值点处的横坐标%f 求得的拉格朗日插值多项式的值 n=length(X);newt=[X',Y'];%计算差商表for j=2:nfor i=n:-1:1if i>=jY(i)=(Y(i)-Y(i-1))/(X(i)-X(i-j+1)); else Y(i)=0;endendnewt=[newt,Y'];end%计算牛顿插值f=newt(1,2);for i=2:nz=1;for k=1:i-1z=(xi-X(k))*z;endf=f+newt(i-1,i)*z;endfprintf('%d\n',f)return3三次样条插值第一类边界条件Threch.mfunction S=Threch1(X,Y,dy0,dyn,xi)% X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%dy0左端点处的一阶导数% dyn右端点处的一阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求函数的一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求函数的二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1);%¸赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))...+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));%三样条插值函数表达式endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M (i+1)-M(i))/(6*h(i))*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return4 三次样条插值第二类边界条件Threch2.mfunction [Sx]=Threch2(X,Y,d2y0,d2yn,xi)X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%d2y0左端点处的二阶导数% d2yn右端点处的二阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=2*d2y0;d(n+1)=2*d2yn;%赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=0;A(n+1,n)=0;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2 +(M(i+1)-M(i))/(6*h(i))*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return5插值节点处的插值结果main3.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];xi=0.13;%xi=0.36;disp('xi=0.13');%disp('xi=0.36');disp('拉格朗日插值结果');lang(X,Y,xi);disp('牛顿插值结果');newton(X,Y,xi);disp('三次样条第一类边界条件插值结果');Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分别为两端点处的一阶导数disp('三次样条第二类边界条件插值结果');Threch2(X,Y,0,-0.136,xi);%0,-0.136分别为两端点处的二阶导数6将多种插值函数即原函数图像画在同一图上main2.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];a=linspace(0,0.4,21);NUM=21;L=zeros(1,NUM);N=zeros(1,NUM);S=zeros(1,NUM);B=zeros(1,NUM);for i=1:NUMxi=a(i);L(i)=lang(X,Y,xi);% 拉格朗日插值N(i)=newton(X,Y,xi);%牛顿插值B(i)=normcdf(xi,0,1);%原函数S(i)=Threch1(X,Y,0.4,0.36,xi);%三次样条函数第一类边界条件endplot(a,B,'--r');hold on;plot(a,L,'b');hold on;plot(a,N,'r');hold on;plot(a,S,'r+');hold on;legend('原函数','拉格朗日插值','牛顿插值','三次样条插值',2);hold off7增加插值节点观察误差变化main4.mclear;clc;N=5;%4.5第一问Ini=zeros(1,1001);a=linspace(-1,1,1001);Ini=1./(1+25*a.^2);for i=1:3 %节点数量变化次数N=2*N;t=linspace(-1,1,N+1);%插值节点ft=1./(1+25*t.^2);%插值节点函数值val=linspace(-1,1,101);for j=1:101L(j)=lang(t,ft,val(j));S(j)=Threch1(t,ft,0.074,-0.074,val(j));%三样条第一类边界条件插值endplot(a,Ini,'k')%原函数图象hold onplot(val,L,'r')%拉格朗日插值函数图像hold onplot(val,S,'b')%三次样条插值函数图像str=sprintf('插值节点为%d时的插值效果',N);title(str);legend('原函数','拉格朗日插值','三次样条插值');%显示图例hold offfigureend8车门曲线main5.mclearclcX=[0,1,2,3,4,5,6,7,8,9,10];Y=[0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29]; dy0=0.8;dyn=0.2;n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:nh(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:nf2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1); A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;x=zeros(1,n);S=zeros(1,n);for i=1:nx(i)=X(i)+0.5;S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x(i)-X(i))+M(i)/2*(x(i)-X(i ))^2+(M(i+1)-M(i))/(6*h(i))*(x(i)-X(i))^3;endplot(X,Y,'k'); hold on;plot(x,S,'o');title('三次样条插值效果图');legend('已知插值节点','三次样条插值');hold off实验结果：4.31计算插值节点处的函数值xi=0.13时Xi=0.36时2将多种插值函数即原函数图像画在同一图上4.5.1增加插值节点观察误差变化从上面三图可以看出增加插值节点并不能改善差之效果4.5.2 车门曲线。
MATLAB大作业 给定一个时间序列,使用三次样条插值方法进行均匀内插
MATLAB作业给定一个时间序列，使用三次样条插值方法进行均匀内插（题目的相关说明：按题目要求编写一个MATLAB程序函数，并把自己编制程序所得的结果与MATLAB库函数分析结果进行对比。
）理论基础：时间序列的概念：时间序列是一种定量预测方法，又称简单外延法，时间序列分析是根据系统观测得到的时间序列数据，通过曲线拟合和参数估计来建立数学模型的理论与方法，时间序列分析可分为以下三种情况（1）把一个时间序列的值变动为N 个组成部分，通常可以分为四种 a、倾向变动，又称长期趋势变动 b、循环变动，又称周期变动 c、季节变动，即每年有规则的反复进行变动 d、不规则变动，即随机变动。
然后把这四个综合到一起得出预测的结果。
虽然分成这四部分，但这四部分之间的相互关系是怎么样的呢，目前一般采用相乘的关系，其实各个部分都是在其他部分作用的基础上进行作用的，所以采用相乘是有一定依据的，此种方法适合于短期预测和库存预测（2）把预测对象、预测目标和对预测的影响因素都看成为具有时序的，为时间的函数，而时间序列法就是研究预测对象自身变化过程及发展趋势，如果未来预测是线性的，其数学模型为YT+L=aT+bTL,YT+L为未来预测值，aT为截距，bT为斜率，L为由T到需要预测的单位时间数（如5年、10年等）（3）根据预测对象与影响因素之间的关系及影响程度来推算未来，与目标的相关因素很多，只能选择那些因果关系较强的为预测影响的因素，此时间序列法用于短期预测比较有效，若要用于长期预测，还需要结合其他方法才行。
三次样条插值的实际应用：在制造船体和汽车外形等工艺中传统的设计方法是，首先由设计人员按外形要求，给出外形曲线的一组离散点值，施工人员准备好有弹性的样条（一般用竹条或有弹性的钢条）和压铁，将压铁放在点的位置上，调整竹条的形状，使其自然光滑，这时竹条表示一条插值曲线，我们称为样条函数。
从数学上看，这一条近似于分段的三次多项式，在节点处具有一阶和二阶连续微商。
计算方法上机作业
计算方法上机报告上机实习题目1.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测(1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；（1）算法思想分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。
在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。
计算光缆长度时，用如下公式：20()L f x ds =⎰200(f x =⎰ 191(k kk f x +==∑⎰= 本题采取三次样条插值的方法，因为三次样条插值方法是一种具有较好“光滑性”的分段插值方法。
根据提供的数据，只用x,y 值，不包含导数值，因此采用第三类三次插值多项式进行插值编程。
设计算法如下：1. For n i ,,2,1,0⋅⋅⋅= 1.1 i i M y ⇒2. For 2,1=k2.1 For k n n i ,,1, -=2.1.1 i k i i i i M x x M M ⇒----)/()(13. 101h x x ⇒-4. For 1-,,2,1n i = 4.1 11++⇒-i i i h x x4.2 b a c c h h h i i i i i i ⇒⇒-⇒+++2;1;)/(11 4.3 i i d M ⇒+165. 0000;;c M d M d n n ⇒⇒⇒λn n n b a b ⇒⇒⇒2;;20μ6. 1111,γμ⇒⇒d b7. 获取M 的矩阵元素个数，存入m8. For m k ,,3,2 = 8.1 k k k l a ⇒-1/μ 8.2 k k k k c l b μ⇒⋅-1- 8.3 k k k k l d γγ⇒⋅-1- 9. m m m M ⇒μγ/10. For 1,,2,1 --=m m k 10.1 k k k k k M M c ⇒⋅-+μγ/)(1 11. 获取x 的元素个数存入s 12.k ⇒113. For 1,,2,1-=s i13.1 if i x x ≤~then k i ⇒；break else k i ⇒+114.xx x x x x h x x k k k k ˆ~;~;11⇒-⇒-⇒--- y h x h M y x h M y x M x M k k k k k k ~/]ˆ)6()6(6ˆ6[2211331⇒-+-++---（3）源程序clear; clc;x=0:1:20; %产生从0到20含21个等分点的数组 X=0:0.2:20;y=[9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93]; %等分点位置的深度数据 n=length(x); %等分点的数目N=length(X);%% 求三次样条插值函数s(x)M=y;for k=2:3; %计算二阶差商并存放在M中 for i=n:-1:k;M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));endendh(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量d for i=2:n-1;h(i)=x(i+1)-x(i);c(i)=h(i)/(h(i)+h(i-1));a(i)=1-c(i);b(i)=2;d(i)=6*M(i+1);endM(1)=0; %选择自然边界条件M(n)=0;b(1)=2;b(n)=2;c(1)=0;a(n)=0;d(1)=0;d(n)=0;u(1)=b(1); %对三对角阵进行LU分解y1(1)=d(1);for k=2:n;l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*c(k-1);y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n); %追赶法求解样条参数M(i)for k=n-1:-1:1;M(k)=(y1(k)-c(k)*M(k+1))/u(k);ends=zeros(1,N);for m=1:N;k=1;for i=2:n-1if X(m)<=x(i);k=i-1;break;elsek=i;endendH=x(k+1)-x(k); %在各区间用三次样条插值函数计算X点处的值x1=x(k+1)-X(m);x2=X(m)-x(k);s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-(M(k)*(H^2)/6)) *x1+(y(k+1)-(M(k+1)*(H^2)/6))*x2)/H;end%% 计算所需光缆长度L=0; %计算所需光缆长度for i=2:NL=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);enddisp('所需光缆长度为 L=');disp(L);figureplot(x,y,'*',X,s,'-') %绘制铺设河底光缆的曲线图xlabel('位置X/测量点','fontsize',16); %标注坐标轴含义ylabel('深度Y/m','fontsize',16);title('铺设河底光缆的曲线图','fontsize',16);grid;（4）结果与分析铺设海底光缆的曲线图如下图所示：拟合结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图。
matlab_牛顿插值法_三次样条插值法
(){}21()(11),5,10,20:12521()1,(0,1,2,,)()2,(0,1,2,,)()()235,20:1100(i i ii n n k k k Newton f x x n x f x x i i n f x nxy i n Newton N x S x n x k y f x =-≤≤=+=-+====-+= 题目：插值多项式和三次样条插值多项式。
已知对作、计算函数在点处的值；、求插值数据点的插值多项式和三次样条插值多项式；、对计算和相应的函数值),()() (1,2,,99)4:()max ()()max()n k n k n k n k n k n k kkN x S x k E N y N x E S y S x ==-=- 和；、计算，；解释你所得到的结果。
算法组织:本题在算法上需要解决的问题主要是：求出第二问中的Newton 插值多项式)(x N n 和三次样条插值多项式()n S x 。
如此，则第三、四问则迎刃而解。
计算两种插值多项式的算法如下：一、求Newton 插值多项式)(x N n ，算法组织如下：Newton 插值多项式的表达式如下：)())(()()(110010--⋅⋅⋅--+⋅⋅⋅+-+=n n n x x x x x x c x x c c x N其中每一项的系数c i 的表达式如下：1102110),,,(),,,(),,,(x x x x x f x x x f x x x f c i i i i i -⋅⋅⋅-⋅⋅⋅=⋅⋅⋅=-根据i c 以上公式，计算的步骤如下：⎪⎪⎪⎩⎪⎪⎪⎨⎧⋅⋅⋅+⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅----),,,,(1),,,(),,,,(),(,),,(2)(,),(),(11101111011010n n n n n n n n x x x x f n x x x f x x x f n x x f x x f x f x f x f 、计算、计算、计算、计算 二、求三次样条插值多项式)(x S n ，算法组织如下：所谓三次样条插值多项式)(x S n 是一种分段函数，它在节点i x 011()n n a x x x x b -=<<⋅⋅⋅<<=分成的每个小区间1[,]i i x x -上是3次多项式，其在此区间上的表达式如下：22331111111()[()()]()()666[,]1,2,,.i i i i i i i i i i i i i i i i i h x x h x x S x x x M x x M y M y M h h h x x x i n --------=-+-+-+-∈=⋅⋅⋅，，因此，只要确定了i M 的值，就确定了整个表达式，i M 的计算方法如下： 令：11111111116()6(,,)i i i i i i i i i i i i i ii i i i i i i h h h h h h y y y y d f x x x h h h h μλμ++++--+++⎧===-⎪++⎪⎨--⎪=-=⎪+⎩， 则i M 满足如下n-1个方程：1121,2,,1i i i i i i M M M d i n μλ-+++==⋅⋅⋅-，方程中有n+1个未知量，则令0M 和n M 分别为零，则由上面的方程组可得到(11)i M i n ≤≤-的值，可得到整个区间上的三次样条插值多项式)(x S n 。
(完整word版)数值分析作业-三次样条插值..
数值计算方法作业实验4.3 三次样条差值函数实验目的：掌握三次样条插值函数的三弯矩方法。
实验函数:dt ex f xt ⎰∞--=2221)(πx 0.0 0.1 0.2 0.3 0.4 F(x) 0.5000 0.5398 0.57930.61790.7554求f(0.13)和f(0.36)的近似值实验内容：（1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值；（3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线比较插值结果。
实验4.5 三次样条差值函数的收敛性实验目的：多项式插值不一定是收敛的，即插值的节点多，效果不一定好。
对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。
实验内容：按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。
实验要求：（1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；（2） 三次样条插值函数的思想最早产生于工业部门。
作为工业应用的例子，考实验名称 实验 4.3三次样条插值函数（P126）4.5三次样条插值函数的收敛性（P127） 实验时间姓名班级学号成绩虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一段数据如下： k x 0 1 2 3 4 5 6 7 8 9 10 k y 0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 ky ' 0.80.2算法描述：拉格朗日插值：错误！未找到引用源。
其中错误！未找到引用源。
是拉格朗日基函数，其表达式为：()∏≠=--=ni j j j i ji x x x x x l 0)()(牛顿插值：))...()(](,...,,[....))(0](,,[)0](,[)()(1102101210100----++--+-+=n n n x x x x x x x x x x f x x x x x x x f x x x x f x f x N其中⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧--=--=--=-)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[01102110x x x x x f x x x f x x x f x x x x f x x f x x x f x x x f x f x x f n n n n i k j i k j k j i ji j i j i三样条插值：所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a<X0<X1……<Xn<b)分成的每个小区间[x i-1,x i ]上是三次多项式，其在此区间上的表达式如下：],[),6()6(]6)([6)(6)()(111113131i i ii i i i i i i i i i i i i i i i i i x x x h yM h M h h y x M M h h y y h x x Mi h x x M x S -------∈-+-+---+-+-=式中Mi=)(i x S ''.因此，只要确定了Mi 的值，就确定了整个表达式，Mi 的计算方法如下：令⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=---+=+=+=+--++++++],,[6)(6111111111i i i i i i i i i i i i i i i i i i i ix x x f h y y h y y h h d h h h h h h λμ则Mi 满足如下n-1个方程：1,...2,1,211-==+++-n i d M M M i i i i i i λμ 常用的边界条件有如下几类：（1） 给定区间两端点的斜率m 0,m n ,即n n n m y x S m y x S ='='='=')(,)(000 （2） 给定区间两端点的二阶导数M0，Mn,即n n n M y x S M y x S =''=''=''='')(,)(000 （3） 假设y=f(x)是以b-a 为周期的周期函数，则要求三次样条插值函数S （x ）也为周期函数，对S （x ）加上周期条件2,1,0),0()0()(0)(=-=+p x S x S n p p对于第一类边界条件有⎪⎪⎩⎪⎪⎨⎧--=+--=+--)(62)(6211001110n n n n n n i h y y mn h M M m h y y h M M对于第二类边界条件有⎩⎨⎧=+=+-n n n n d M M d M M 221100μλ其中n n n n nnn M u x x f m h d M m x x f h d )1(2]),[(6)1(2)],[(6100001010-+-=-+-=-μλλ那么解就可以为⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n n n n d d d d d M M M M M 1210121011...2...............2............................1..2.1......0..2μλμλμλ 对于第三类边界条件，)0()0(,,000-=+==n n n x S x S M M y y ，由此推得0010012d M M M n =-++μλ，其中]),1[],[(6,,101010110n n nn n n x x f x x f h h d h h h h h h --+=+=+=μλ，那么解就可以为： ⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-------1221012101221100...2.............2..............................2..,,.......,..22n n n n n n n d d d d d M M M M M n μλλμλμμλ 程序代码： 1拉格朗日插值函数Lang.mfunction f=lang(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标 %xi 插值点处的横坐标%f 求得的拉格朗日插值多项式的值 n=length(X); f=0; for i=1:n l=1; for j=1:i-1l=l.*(xi-X(j))/(X(i)-X(j)); end ; for j=i+1:nl=l.*(xi-X(j))/(X(i)-X(j)); end ;%拉格朗日基函数 f=f+l*Y(i); endfprintf('%d\n',f) return2 牛顿插值函数newton.mfunction f=newton(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标 %xi 插值点处的横坐标%f 求得的拉格朗日插值多项式的值 n=length(X);newt=[X',Y'];%计算差商表for j=2:nfor i=n:-1:1if i>=jY(i)=(Y(i)-Y(i-1))/(X(i)-X(i-j+1));else Y(i)=0;endendnewt=[newt,Y'];end%计算牛顿插值f=newt(1,2);for i=2:nz=1;for k=1:i-1z=(xi-X(k))*z;endf=f+newt(i-1,i)*z;endfprintf('%d\n',f)return3三次样条插值第一类边界条件Threch.mfunction S=Threch1(X,Y,dy0,dyn,xi)% X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%dy0左端点处的一阶导数% dyn右端点处的一阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求函数的一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求函数的二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1);%¸赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))...+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));%三样条插值函数表达式endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M (i+1)-M(i))/(6*h(i))*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return4 三次样条插值第二类边界条件Threch2.mfunction [Sx]=Threch2(X,Y,d2y0,d2yn,xi)X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%d2y0左端点处的二阶导数% d2yn右端点处的二阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=2*d2y0;d(n+1)=2*d2yn;%赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=0;A(n+1,n)=0;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2 +(M(i+1)-M(i))/(6*h(i))*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return5插值节点处的插值结果main3.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];xi=0.13;%xi=0.36;disp('xi=0.13');%disp('xi=0.36');disp('拉格朗日插值结果');lang(X,Y,xi);disp('牛顿插值结果');newton(X,Y,xi);disp('三次样条第一类边界条件插值结果');Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分别为两端点处的一阶导数disp('三次样条第二类边界条件插值结果');Threch2(X,Y,0,-0.136,xi);%0,-0.136分别为两端点处的二阶导数6将多种插值函数即原函数图像画在同一张图上main2.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];a=linspace(0,0.4,21);NUM=21;L=zeros(1,NUM);N=zeros(1,NUM);S=zeros(1,NUM);B=zeros(1,NUM);for i=1:NUMxi=a(i);L(i)=lang(X,Y,xi);% 拉格朗日插值N(i)=newton(X,Y,xi);%牛顿插值B(i)=normcdf(xi,0,1);%原函数S(i)=Threch1(X,Y,0.4,0.36,xi);%三次样条函数第一类边界条件endplot(a,B,'--r');hold on;plot(a,L,'b');hold on;plot(a,N,'r');hold on;plot(a,S,'r+');hold on;legend('原函数','拉格朗日插值','牛顿插值','三次样条插值',2);hold off7增加插值节点观察误差变化main4.mclear;clc;N=5;%4.5第一问Ini=zeros(1,1001);a=linspace(-1,1,1001);Ini=1./(1+25*a.^2);for i=1:3 %节点数量变化次数N=2*N;t=linspace(-1,1,N+1);%插值节点ft=1./(1+25*t.^2);%插值节点函数值val=linspace(-1,1,101);for j=1:101L(j)=lang(t,ft,val(j));S(j)=Threch1(t,ft,0.074,-0.074,val(j));%三样条第一类边界条件插值endplot(a,Ini,'k')%原函数图象hold onplot(val,L,'r')%拉格朗日插值函数图像hold onplot(val,S,'b')%三次样条插值函数图像str=sprintf('插值节点为%d时的插值效果',N);title(str);legend('原函数','拉格朗日插值','三次样条插值');%显示图例hold offfigureend8车门曲线main5.mclearclcX=[0,1,2,3,4,5,6,7,8,9,10];Y=[0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29]; dy0=0.8;dyn=0.2;n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:nh(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:nf2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1); A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;x=zeros(1,n);S=zeros(1,n);for i=1:nx(i)=X(i)+0.5;S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x(i)-X(i))+M(i)/2*(x(i)-X(i ))^2+(M(i+1)-M(i))/(6*h(i))*(x(i)-X(i))^3;endplot(X,Y,'k'); hold on;plot(x,S,'o');title('三次样条插值效果图');legend('已知插值节点','三次样条插值');hold off实验结果：4.31计算插值节点处的函数值xi=0.13时Xi=0.36时2将多种插值函数即原函数图像画在同一张图上4.5.1增加插值节点观察误差变化4.5.2 车门曲线。
数值分析作业-三次样条插值
数值计算方法作业实验4.3 三次样条差值函数实验目的：掌握三次样条插值函数的三弯矩方法。
实验函数:dt ex f xt ⎰∞--=2221)(π实验内容：（1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值；（3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线比较插值结果。
实验4.5 三次样条差值函数的收敛性实验目的：多项式插值不一定是收敛的，即插值的节点多，效果不一定好。
对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。
实验内容：按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。
实验要求：（1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；（2） 三次样条插值函数的思想最早产生于工业部门。
作为工业应用的例子，考虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一段数据如下：kx012345678910 ky0.00.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29ky'0.80.2算法描述：拉格朗日插值：其中是拉格朗日基函数，其表达式为：()∏≠=--=nijj jiji xxxxxl)()(牛顿插值：))...()(](,...,,[....))(](,,[)0](,[)()(11211211----++--+-+=nnnxxxxxxxxxxfxxxxxxxfxxxxfxfxN其中⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧--=--=--=-)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[11211xxxxxfxxxfxxxfxxxxfxxfxxxfxxxfxfxxfnnnnikjikjkjijijiji三样条插值：所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a<X0<X1……<Xn<b)分成的每个小区间[xi-1,xi]上是三次多项式，其在此区间上的表达式如下：],[),6()6(]6)([6)(6)()(111113131iiiiiiiiiiiiiiiiiiiiixxxhyMhMhhyxMMhhyyhxxMihxxMxS-------∈-+-+---+-+-=式中Mi=)(ixS''.因此，只要确定了Mi 的值，就确定了整个表达式，Mi 的计算方法如下：令⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=---+=+=+=+--++++++],,[6)(6111111111i i i i i i i i i i i i i i i i i i i ix x x f h y y h y y h h d h h h h h h λμ则Mi 满足如下n-1个方程：1,...2,1,211-==+++-n i d M M M i i i i i i λμ 常用的边界条件有如下几类：（1） 给定区间两端点的斜率m 0,m n ,即n n n m y x S m y x S ='='='=')(,)(000 （2） 给定区间两端点的二阶导数M0，Mn,即n n n M y x S M y x S =''=''=''='')(,)(000 （3） 假设y=f(x)是以b-a 为周期的周期函数，则要求三次样条插值函数S （x ）也为周期函数，对S （x ）加上周期条件2,1,0),0()0()(0)(=-=+p x S x S n p p对于第一类边界条件有⎪⎪⎩⎪⎪⎨⎧--=+--=+--)(62)(6211001110n n n n n n i h y y mn h M M m h y y h M M对于第二类边界条件有⎩⎨⎧=+=+-n n n n d M M d M M 221100μλ其中n n n n nnn M u x x f m h d M m x x f h d )1(2]),[(6)1(2)],[(6100001010-+-=-+-=-μλλ那么解就可以为⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n n n n d d d d d M M M M M 1210121011...2...............2............................1..2.1......0..2μλμλμλ对于第三类边界条件，)0()0(,,000-=+==n n n x S x S M M y y ，由此推得0010012d M M M n =-++μλ，其中]),1[],[(6,,101010110n n nn n n x x f x x f h h d h h h h h h --+=+=+=μλ，那么解就可以为： ⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-------1221012101221100...2.............2..............................2..,,.......,..22n n n n n n n d d d d d M M M M M n μλλμλμμλ 程序代码： 1拉格朗日插值函数Lang.mfunction f=lang(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标 %xi 插值点处的横坐标%f 求得的拉格朗日插值多项式的值 n=length(X); f=0; for i=1:n l=1; for j=1:i-1l=l.*(xi-X(j))/(X(i)-X(j)); end ; for j=i+1:nl=l.*(xi-X(j))/(X(i)-X(j)); end ;%拉格朗日基函数 f=f+l*Y(i); endfprintf('%d\n',f) return2 牛顿插值函数newton.mfunction f=newton(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标%xi插值点处的横坐标%f求得的拉格朗日插值多项式的值n=length(X);newt=[X',Y'];%计算差商表for j=2:nfor i=n:-1:1if i>=jY(i)=(Y(i)-Y(i-1))/(X(i)-X(i-j+1));else Y(i)=0;endendnewt=[newt,Y'];end%计算牛顿插值f=newt(1,2);for i=2:nz=1;for k=1:i-1z=(xi-X(k))*z;endf=f+newt(i-1,i)*z;endfprintf('%d\n',f)return3三次样条插值第一类边界条件Threch.mfunction S=Threch1(X,Y,dy0,dyn,xi)% X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%dy0左端点处的一阶导数% dyn右端点处的一阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求函数的一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求函数的二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1);%¸赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))...+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));%三样条插值函数表达式endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6 *h(i))*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return4 三次样条插值第二类边界条件Threch2.mfunction [Sx]=Threch2(X,Y,d2y0,d2yn,xi)X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%d2y0左端点处的二阶导数% d2yn右端点处的二阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=2*d2y0;d(n+1)=2*d2yn;%赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=0;A(n+1,n)=0;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i)) /(6*h(i))*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return5插值节点处的插值结果main3.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];xi=0.13;%xi=0.36;disp('xi=0.13');%disp('xi=0.36');disp('拉格朗日插值结果');lang(X,Y,xi);disp('牛顿插值结果');newton(X,Y,xi);disp('三次样条第一类边界条件插值结果');Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分别为两端点处的一阶导数disp('三次样条第二类边界条件插值结果');Threch2(X,Y,0,-0.136,xi);%0,-0.136分别为两端点处的二阶导数6将多种插值函数即原函数图像画在同一张图上main2.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];a=linspace(0,0.4,21);NUM=21;L=zeros(1,NUM);N=zeros(1,NUM);S=zeros(1,NUM);B=zeros(1,NUM);for i=1:NUMxi=a(i);L(i)=lang(X,Y,xi);% 拉格朗日插值N(i)=newton(X,Y,xi);%牛顿插值B(i)=normcdf(xi,0,1);%原函数S(i)=Threch1(X,Y,0.4,0.36,xi);%三次样条函数第一类边界条件endplot(a,B,'--r');hold on;plot(a,L,'b');hold on;plot(a,N,'r');hold on;plot(a,S,'r+');hold on;legend('原函数','拉格朗日插值','牛顿插值','三次样条插值',2);hold off7增加插值节点观察误差变化main4.mclear;clc;N=5;%4.5第一问Ini=zeros(1,1001);a=linspace(-1,1,1001);Ini=1./(1+25*a.^2);for i=1:3 %节点数量变化次数N=2*N;t=linspace(-1,1,N+1);%插值节点ft=1./(1+25*t.^2);%插值节点函数值val=linspace(-1,1,101);for j=1:101L(j)=lang(t,ft,val(j));S(j)=Threch1(t,ft,0.074,-0.074,val(j));%三样条第一类边界条件插值endplot(a,Ini,'k')%原函数图象hold onplot(val,L,'r')%拉格朗日插值函数图像hold onplot(val,S,'b')%三次样条插值函数图像str=sprintf('插值节点为%d时的插值效果',N);title(str);legend('原函数','拉格朗日插值','三次样条插值');%显示图例hold offfigureend8车门曲线main5.mclearclcX=[0,1,2,3,4,5,6,7,8,9,10];Y=[0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29]; dy0=0.8;dyn=0.2;n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:nh(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:nf2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1); A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;x=zeros(1,n);S=zeros(1,n);for i=1:nx(i)=X(i)+0.5;S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x(i)-X(i))+M(i)/2*(x(i)-X(i))^2+(M(i+1)-M (i))/(6*h(i))*(x(i)-X(i))^3;endplot(X,Y,'k'); hold on;plot(x,S,'o');title('三次样条插值效果图');legend('已知插值节点','三次样条插值');hold off实验结果：4.31计算插值节点处的函数值xi=0.13时Xi=0.36时2将多种插值函数即原函数图像画在同一张图上4.5.1增加插值节点观察误差变化从上面三张图可以看出增加插值节点并不能改善差之效果4.5.2 车门曲线（注：专业文档是经验性极强的领域，无法思考和涵盖全面，素材和资料部分来自网络，供参考。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
例1 设)(x f 为定义在[0,3]上的函数，有下列函数值表：且2.0)('0=x f ，1)('3-=x f ，试求区间[0,3]上满足上述条件的三次样条插值函数)(x s本算法求解出的三次样条插值函数将写成三弯矩方程的形式：)()6()()6()(6)(6)(211123131j j jj j j jj j j j jj j jj x x h h M y x x h h M y x x h M x x h M x s --+--+-+-=+++++其中，方程中的系数jj h M 6，jj h M 61+，jj j j h h M y )6(2-，jjj j h h M y )6(211++-将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。
以下为Matlab 代码：%=============================% 本段代码解决作业题的例1%============================= clear all clc% 自变量x 与因变量y ，两个边界条件的取值 IndVar = [0, 1, 2, 3]; DepVar = [0, 0.5, 2, 1.5];LeftBoun = 0.2; RightBoun = -1;% 区间长度向量，其各元素为自变量各段的长度 h = zeros(1, length(IndVar) - 1); for i = 1 : length(IndVar) - 1h(i) = IndVar(i + 1) - IndVar(i); end% 为向量μ赋值mu = zeros(1, length(h));for i = 1 : length(mu) - 1mu(i) = h(i) / (h(i) + h(i + 1));endmu(i + 1) = 1;% 为向量λ赋值lambda = zeros(1, length(h));lambda(1) = 1;for i = 2 : length(lambda)lambda(i) = h(i) / (h(i - 1) + h(i));end% 为向量d赋值d = zeros(1, length(h) + 1);d(1) = 6 * ( (DepV ar(2) - DepVar(1) ) / ( IndVar(2) - IndVar(1) ) - LeftBoun) / h(1);for i = 2 : length(h)a = ( DepVar(i) - DepVar(i - 1) ) / ( IndVar(i) - IndVar(i - 1) );b = ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) );c = (b - a) / ( IndVar(i + 1) - IndVar(i - 1) );d(i) = 6 * c;endd(i + 1) = 6 *( RightBoun - ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) ) ) / h(i);% 为矩阵A赋值% 将主对角线上的元素全部置为2A = zeros( length(d), length(d) );for i = 1 : length(d)A(i, i) = 2;end% 将向量λ的各元素赋给主对角线右侧第一条对角线for i = 1 : length(d) - 1A(i, i + 1) = lambda(i);end% 将向量d的各元素赋给主对角线左侧第一条对角线for i = 1 : length(d) - 1A(i + 1, i) = mu(i);end% 求解向量MM =A \ d';% 求解每一段曲线的函数表达式for i = 1 : length(h)Coefs_1 = M(i) / (6 * h(i));Part_1 = conv( Coefs_1, ...conv( [-1, IndVar(i + 1)], ...conv( [-1, IndVar(i + 1)], [-1, IndVar(i + 1)] ) ) );S_1 = polyval (Part_1, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_2 = M(i + 1)/(6 * h(i));Part_2 = conv( Coefs_2, ...conv( [1, -IndVar(i)], ...conv( [1, -IndVar(i)], [1, -IndVar(i)] ) ) );S_2 = polyval (Part_2, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_3 = (DepVar(i) - M(i) * h(i)^2 / 6) / h(i);Part_3 = conv(Coefs_3, [-1, IndVar(i + 1)]);S_3 = polyval (Part_3, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_4 = (DepVar(i + 1) - M(i + 1) * h(i)^2 / 6) / h(i);Part_4 = conv(Coefs_4, [1, -IndVar(i)]);S_4 = polyval (Part_4, [IndVar(i) : 0.01 : IndVar(i + 1)]);S = S_1 + S_2 + S_3 + S_4;plot ([IndVar(i) : 0.01 : IndVar(i + 1)], S, 'LineWidth', 1.25)% 在样条插值曲线的相应位置标注该段曲线的函数表达式text(i - 1, polyval(Part_1, 3), ...['\itS', num2str(i), '(x)=', num2str(Coefs_1), '(', num2str( IndVar(i + 1) ), '-x)^{3}+', ...num2str(Coefs_2), '(x-', num2str( IndVar(i) ), ')^{3}+', num2str(Coefs_3), ...'(', num2str( IndVar(i + 1) ), '-x)+', num2str(Coefs_4), '(x-', num2str( IndVar(i) ), ')'], ...'FontName', 'Times New Roman', 'FontSize', 14)hold onend% 过x=1和x=2两个横轴点作垂线%line([1, 1], [2.5, -0.5], 'LineStyle', '--');line([2, 2], [2.5, -0.5], 'LineStyle', '--');% 为x轴和y轴添加标注xlabel( '\itx', 'FontName', 'Times New Roman', ...'FontSize', 14, 'FontWeight', 'bold');ylabel( '\its(x)', 'FontName', 'Times New Roman', ...'Rotation', 0, 'FontSize', 14, 'FontWeight', 'bold');最终，三次样条插值函数s(x)表达式为：[][][]⎪⎩⎪⎨⎧∈-+-+-+--∈-+-+---∈+-++--=.3,2,)2(44.1)3(62.2)2(06.0)3(62.0,2,1,)1(62.2)2(08.0)1(62.0)2(42.0,1,0,08.0)1(06.042.0)1(06.0)(333333x x x x x x x x x x x x x x x x s曲线的图像如图所示：例2 已知函数值表：试求在区间[1,5]上满足上述函数表所给出的插值条件的三次自然样条插值函数)(x s本算法求解出的三次样条插值函数将写成三弯矩方程的形式：)()6()()6()(6)(6)(211123131j j jj j j jj j j j jj j jj x x h h M y x x h h M y x x h M x x h M x s --+--+-+-=+++++其中，方程中的系数jj h M 6，jj h M 61+，jj j j h h M y )6(2-，jjj j h h M y )6(211++-将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。
以下为Matlab 代码：%=============================% 本段代码解决作业题的例2%============================= clear all clc% 自变量x 与因变量y 的取值 IndVar = [1, 2, 4, 5]; DepVar = [1, 3, 4, 2];% 区间长度向量，其各元素为自变量各段的长度 h = zeros(1, length(IndVar) - 1); for i = 1 : length(IndVar) - 1h(i) = IndVar(i + 1) - IndVar(i); end% 为向量μ赋值mu = zeros(1, length(h)); for i = 1 : length(mu) - 1mu(i) = h(i) / (h(i) + h(i + 1));endmu(i + 1) = 0;% 为向量λ赋值lambda = zeros(1, length(h));lambda(1) = 0;for i = 2 : length(lambda)lambda(i) = h(i) / (h(i - 1) + h(i));end% 为向量d赋值d = zeros(1, length(h) + 1);d(1) = 0;for i = 2 : length(h)a = ( DepVar(i) - DepVar(i - 1) ) / ( IndVar(i) - IndVar(i - 1) );b = ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) );c = (b - a) / ( IndVar(i + 1) - IndVar(i - 1) );d(i) = 6 * c;endd(i + 1) = 0;% 为矩阵A赋值% 将主对角线上的元素全部置为2A = zeros( length(d), length(d) );for i = 1 : length(d)A(i, i) = 2;end% 将向量λ的各元素赋给主对角线右侧第一条对角线for i = 1 : length(d) - 1A(i, i + 1) = lambda(i);end% 将向量d的各元素赋给主对角线左侧第一条对角线for i = 1 : length(d) - 1A(i + 1, i) = mu(i);end% 求解向量MM =A \ d';% 求解每一段曲线的函数表达式for i = 1 : length(h)Coefs_1 = M(i) / (6 * h(i));Part_1 = conv( Coefs_1, ...conv( [-1, IndVar(i + 1)], ...conv( [-1, IndVar(i + 1)], [-1, IndVar(i + 1)] ) ) );S_1 = polyval (Part_1, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_2 = M(i + 1)/(6 * h(i));Part_2 = conv( Coefs_2, ...conv( [1, -IndVar(i)], ...conv( [1, -IndVar(i)], [1, -IndVar(i)] ) ) );S_2 = polyval (Part_2, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_3 = (DepVar(i) - M(i) * h(i)^2 / 6) / h(i);Part_3 = conv(Coefs_3, [-1, IndVar(i + 1)]);S_3 = polyval (Part_3, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_4 = (DepVar(i + 1) - M(i + 1) * h(i)^2 / 6) / h(i);Part_4 = conv(Coefs_4, [1, -IndVar(i)]);S_4 = polyval (Part_4, [IndVar(i) : 0.01 : IndVar(i + 1)]);S = S_1 + S_2 + S_3 + S_4;plot ([IndVar(i) : 0.01 : IndVar(i + 1)], S, 'LineWidth', 1.25)% 在样条插值曲线的相应位置标注该段曲线的函数表达式text(i, polyval(Part_1, 5), ...['\itS', num2str(i), '(x)=', num2str(Coefs_1), '(', num2str( IndVar(i + 1) ), '-x)^{3}+', ...num2str(Coefs_2), '(x-', num2str( IndVar(i) ), ')^{3}+', num2str(Coefs_3), ...'(', num2str( IndVar(i + 1) ), '-x)+', num2str(Coefs_4), '(x-', num2str( IndVar(i) ), ')'], ...'FontName', 'Times New Roman', 'FontSize', 14)hold onend% 过x=2和x=4两个横轴点作垂线%line([2, 2], [4.5, 0.5], 'LineStyle', '--');line([4, 4], [4.5, 0.5], 'LineStyle', '--');% 为x轴和y轴添加标注xlabel( '\itx', 'FontName', 'Times New Roman', ...'FontSize', 14, 'FontWeight', 'bold');ylabel( '\its(x)', 'FontName', 'Times New Roman', ...'Rotation', 0, 'FontSize', 14, 'FontWeight', 'bold');最终，三次自然样条插值函数s(x)表达式为：[][][]⎪⎩⎪⎨⎧∈-+-+--∈-+-+----∈-+-+--=.5,4,)4(2)5(375.4)5(375.0,4,2,)2(75.2)4(75.1)2(1875.0)4(0625.0,2,1,)1(125.3)2()1(125.0)(3333x x x x x x x x x x x x x x s曲线的图像如图所示：[][][][]⎪⎪⎩⎪⎪⎨⎧∈-+-+--∈-+-+----∈-+-+----∈-+-+--=.53.0,45.0,)45.0(1.9)53.0(3987.8)53.0(1442.2,45.0,39.0,)39.0(1903.11)45.0(417.10)39.0(859.2)45.0(399.2,39.0,30.0,)3.0(9518.6)39.0(1137.6)3.0(5993.1)39.0(4806.3,30.0,25.0,)25.0(9697.10)3.0(10)25.0(2652.6)(333333x x x x x x x x x x x x x x x x x x x s试求在区间[0.25,0.53]上满足上述函数表所给出的插值条件的三次自然样条插值函数s(x)求解出的三次样条插值函数将写成三弯矩方程的形式：)()6()()6()(6)(6)(211123131j jjj j j jj j j j jj j jj x x h h M y x x h h M y x x h M x x h M x s --+--+-+-=+++++本题采用和例2基本相同的Matlab 代码，只改变初始条件。