三次样条插值课后题集
- 格式:doc
- 大小:59.71 KB
- 文档页数:17
问题
分段低次插值
在处理实际问题时,总是希望将所得到的数据点用得越多越好。
最简单的方法是用直线将函数值点直接连接。
分段低次插值
基本思想:用分段低次多项式来代替单个多项式。
具体作法:(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三次样条插值1. 前言在数学和工程领域中,插值是一种常见的数值分析技术,它可以用来估计不连续数据点之间的值。
而三次样条插值作为一种常用的插值方法,在Matlab中有着广泛的应用。
本文将从简单到复杂,由浅入深地解析Matlab中的三次样条插值方法,以便读者更深入地理解这一技术。
2. 三次样条插值概述三次样条插值是一种利用分段三次多项式对数据点进行插值的方法。
在Matlab中,可以使用spline函数来进行三次样条插值。
该函数需要输入数据点的x和y坐标,然后可以根据需要进行插值操作。
3. 三次样条插值的基本原理在进行三次样条插值时,首先需要对数据点进行分段处理,然后在每个分段上构造出一个三次多项式函数。
这些多项式函数需要满足一定的插值条件,如在数据点处函数值相等、一阶导数相等等。
通过这些条件,可以得到一个关于数据点的插值函数。
4. Matlab中的三次样条插值实现在Matlab中,可以使用spline函数来进行三次样条插值。
通过传入数据点的x和y坐标,可以得到一个关于x的插值函数。
spline函数也支持在已知插值函数上进行插值点的求值,这为用户提供了极大的灵活性。
5. 三次样条插值的适用范围和局限性虽然三次样条插值在许多情况下都能够得到较好的插值效果,但也存在一些局限性。
在数据点分布不均匀或有较大噪音的情况下,三次样条插值可能会出现较大的误差。
在实际应用中,需要根据具体情况选择合适的插值方法。
6. 个人观点和总结通过对Matlab中三次样条插值的深度解析,我深刻地理解了这一插值方法的原理和实现方式。
在实际工程应用中,我会根据数据点的情况选择合适的插值方法,以确保得到准确且可靠的结果。
我也意识到插值方法的局限性,这为我在实际工作中的决策提供了重要的参考。
通过以上深度解析,相信读者已经对Matlab中的三次样条插值有了更加全面、深刻和灵活的理解。
在实际应用中,希望读者能够根据具体情况选择合适的插值方法,以提高工作效率和准确性。
例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的值求出。
数值计算方法作业实验4.3三次样条差值函数实验目的:掌握三次样条插值函数的三弯矩方法实验函数:求和的近似值实验内容:(1) 编程实现求三次样条插值函数的算法,分别考虑不同的边界条件;(2) 计算各插值节点的弯矩值;(3) 在同一坐标系中绘制函数f(x),插值多项式,三次样条插值多项式的曲线比较插值结果。
实验4.5三次样条差值函数的收敛性实验目的:多项式插值不一定是收敛的,即插值的节点多,效果不一定好。
对三次样条插值函数如何呢?理论上证明三次样条插值函数的收敛性是比较困难的,通过本实验可以证明这一理论结果。
实验内容:按照一定的规则分别选择等距或非等距的插值节点,并不断增加插值节点的个数。
实验要求:(1)随着节点个数的增加,比较被逼近函数和三样条插值函数的误差变化情况,分析所得结果并与拉格朗日插值多项式比较;(2)三次样条插值函数的思想最早产生于工业部门。
作为工业应用的例子,考虑如下例子:某汽车制造商根据三次样条插值函数设计车门曲线, 其中一 段数据X k 0 1 23 4 5678910y k 0.00.79 1.532.19 2.713.03 3.27 2.89 3.06 3.19 3.29y k0.80.2算法描述:拉格朗日插值:错误!未找到引用源。
n(x _ X ) 其中错误!未找到引用源。
是拉格朗日基函数,其表达式为:h(x)」j=0 (x i- X j )牛顿插值:N n (x) =f (X g ) f[X o ,X i ](X -xO) f[X o ,X i ,X 2〕(X - xO)(x - X i ) •…f[X g ,X i ...X n ] =(f[X i ,X 2,...X n ] - f [ X 。
,为,..人」)/(X . - X g )三样条插值:所谓三次样条插值多项式Sn(x)是一种分段函数,它在节点Xi(a<X0<X1……<Xn<b)分成的每个小区间[x i-i ,x i ]上是三次多项式,其在此区间 上的表达式如下:f[X °,X i ,X 2,...X n ](X -X °)(X -X i )...(X-Xn J )f [X i , X j ]f (X i ) - f (X j ) X i -X jf [X i , X j ,X k]=其中*.f[X j ,X k ] - f[K ,X j ]X k -X iS(x)二 M 3(X i -x) 6h i.Mi (x —Xy )3 . [ y i - y i4 h i (M i - My)6h i h i 6 h ih i M i 4 h i M iy i- 6)( 6*,皿"]因此,只要确定了 Mi 的值,就确定了整个表达式,Mi 的计算方法如下:i 4式中 Mi= S (X i ).则Mi 满足如下n-1个方程:7 M i 」■ 2M i …冷 M i i = di , i =1,2,...n —'1 常用的边界条件有如下几类:(1)给定区间两端点的斜率 m o ,m n ,即s(x 0) = y 0 =m 0,S(x n ) = y n = m n(2) 给定区间两端点的二阶导数 MO ,Mn,即S (XcH y 。
例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 j j 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));% 为向量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 on end% 过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 j j 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 代码,只改变初始条件。