三次样条插值法《数值分析》上机实验作业
昆明理工大学研究生《数值分析》上机实验作业
姓名：
学号：
专业：
一、 课题名称
课题七 三次样条插值法
1、问题提出 设已知数据如下：
求f(x)的三次样条插值函数S(x)。 2、要求
（1）满足自然边界条件0)0.1()2.0(=''=''s s ；
（2）满足第一类边界条件20271.0)2.0(='s ，55741.1)0.1(='s 。 （3）打印输出用追赶法解出的弯曲向量(0M ,1M ,2M ,3M ,4M )和
)1.02.0(i S +（i=0,1,2,3,4,5,6,7,8）的值。并画出)(x S y =的图形。
二、 班级、姓名、学号
三、目的和意义
由于航空、造船等工程设计的需要而发展起来所谓样条插值方法，既保留了分段低次插值多项式的各种优点，又提高了插值函数的光滑性，而且具有较好的稳定性。今天，样条插值方法已成为数值逼近的一个极其重要的分支，在许多领域里得到越来越广泛地应用。其中，尤以三次样条插值函数应用最为广泛，如在高速飞机的机翼形体和船体放样等方面的应用，同时在计算机作图方面更是大有作为。它能够解决一些既有二阶光滑度，又有二阶连续导数的方程，具有良好的收
敛性和稳定性。
1. 通过本次实验进一步了解三次样条插值函数，并通过求解三 弯矩方程组得出曲线函数组；
2. 通过MATLAB 编程实现求三次样条插值函数的算法，分别考虑 不同的边界条件，同时用追赶法解出弯曲向量和)1.02.0(i S +（i=0,1,2,3,4,5,6,7,8）的值。 四、计算公式
首先我们利用)(x S 的二阶导数值),,2,1,0()(j n j M x S j =='
'表达)(x S ,因为在区间]x ,[x 1j j +上)()(j x S x S =是不高于三次的多项式，其二阶导数)(x S ''必是线性函数，所以可表示为：
]x ,[x x ,h x x M h x x M (x)S 1j j j
j
1
j j
1j j
+++∈-+-=''对)(x S ''积分两次并
利用1j 1j j j y )(x S ,y )(x S ++==，可定出积分 常数，于是得三次样条表达式。
j
j
2
j
1j 1j j
1j 2
j j j j
3
j 1
j j
3
1j j h x x )
h 6
M (y h x x )
h 6
M (y 6h )x (x M 6h )x (x M (x)S --
+--
+-+-=+++++
1n ,0,1,2,，j -=
这里),,1,0(n j M j =是未知的。为了确定),,1,0(n j M j =，对)(x S 求导得
)M (M 6
h h y y 2h )x (x M 2h x)(x M (x)S j 1j j j
j
1j j
2
j 1
j j
2
1j j
---+
-+--='++++
由此可得
j j 1j j j
j
1j j M 3
h M 6
h h y y 0)(x S -
-
-=
+'++
类似的可求出)(x S 在区间]x ,[x j 1-j 上的表达式，进而得
j 1j 1j 1j 1
j 1j j j M 3h M 6
h h y y 0)(x S -----+
+-=
-'
利用0)(x S 0)(x S j j +'=-'可得
j 1j j j 1j j d M λ2M M μ=+++- ，1,n ,1,j -= （三弯矩方程）
其中 j
1j 1j j h h h μ+=
--，j
1j j j h h h λ+=
-，
[][
][]
1j j 1
j j
1j j
1j 1j j j x ,x ,x
6f h h x ,x f x ,x f 6
d +---+=+-=.1,n ,1,j -=
其中有（1+n
）个未知数n M M M ,...,10，而方程只有（n-1）个，
当满足第一种边界条件时，可得另两个方程
[]()0100
10,6
2f x x f h M M '-=
+， []()n n n n n n x x f f h M M ,6
211
1----'=
+ 如果令[]()[]()n n n n n n x x f f h g f x x f h d ,6,1,,6
,111
010000---'=='-=
=μλ，将上述方程综合后的一下矩阵形式：
????
???
?????????=??
?
??
?
??
??
?
???
?
????
?????????????----n 1n 10n
1n 10
n 1n 1n 1
10
d d d d M M M M 2μλ2μλ2μλ2 可以证明此方程组满足追赶法的条件，我们用追赶法可得M 的
)1(+n 值，将其带入公式即得)(x s 。
对第二种边界条件，直接的端点方程n n f M f M ''=''=,00并且令
n n n f d f d ''=''===2,2,0000μλ，则又得三弯矩方程同理即可求得解。
五、结构程序设计
1.满足自然边界条件时
自定义函数：followup.m %追赶法求m
%A 为线性方程组的系数矩阵 %b 为常数向量
function m=followup(A,b) n=rank(A); fori=1:n if A(i,i)==0
disp('error:对角元素中有数据为0'); return; end end
d=ones(n,1); a=ones(n-1,1); c=ones(n-1); fori=1:n-1
a(i,1)=A(i+1,i); c(i,1)=A(i,i+1); d(i,1)=A(i,i); end
d(n,1)=A(n,n); fori=2:n
d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1))*c(i-1,1); b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1))*b(i-1,1); end
m(n,1)=b(n,1)/d(n,1); fori=(n-1):-1:1
m(i,1)=(b(i,1)-c(i,1)*m(i+1,1))/d(i,1); end
自定义函数：thrsample2.m %a 为要求的插值点 %f 为区间内的插值函数 %f0为输入点处的插值
%m 为追赶法解出的弯矩向量 function thrsample2(a) x=[0.2:0.2:1.0];
y=[0.9798652 0.9177710 0.8080348 0.6386092 0.3843735];
s02=0;
s10=0;
x0=a;
n=length(x);
fori=1:n
if (x(i)<=x0)&&(x(i+1)>=x0)
index=i;
break;
end
end
A=diag(2*ones(1,n));
A(1,2)=1;
A(n,n-1)=1;
u=zeros(n-2,1);
lamda=zeros(n-1,1);
c=zeros(n,1);
fori=2:n-1
u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));
lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));
c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));
A(i,i+1)=u(i-1);
A(i,i-1)=lamda(i);
end
c(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*s02/2;
c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*s10/2;
m=followup(A,c)
h=x(index+1)-x(index);
syms t;
f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(ind ex))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-m(index+1)*(x(index+1)-t)*(t-x(index)) ^2/h/h
f0=subs(f,'t',x0)
运行结果（方程S、弯矩M和插值函数f的值）为：
S0 = 0.9798652000
S1 = 0.9525726703
S2 = 0.9177710000
S3 = 0.8695049390
S4 = 0.8080348000
S5 = 0.7334085485
S6 = 0.6386092000
S7 = 0.5187304296
S8 =0.3843735000
f =125*((8825841099186633*t)/4503599627370496 - 8825841099186633/45035996273704960)*(t - 2/5)^2 - 125*((8266546267222895*t)/4503599627370496 - 8266546267222895/9007199254740992)*(t - 1/5)^2 -
25*((7396583675403915*t)/18014398509481984 - 1479316735080783/9007199254740992)*(t - 1/5)^2 - 25*(t - 2/5)^2*((2290591133669*t)/8796093022208 - 2290591133669/43980465111040)
2.满足第一类边界条件时
自定义函数：thrsample1.m
%a为要求的插值点
%f为区间内的插值函数
%f0为输入点处的插值
%m为追赶法解出的弯矩向量
function thrsample1(a)
x=[0.2:0.2:1.0];
y=[0.9798652 0.9177710 0.8080348 0.6386092 0.3843735];
s02=0.20271;
s10=1.55741;
x0=a;
n=length(x);
fori=1:n
if (x(i)<=x0)&&(x(i+1)>=x0)
index=i;
break;
end
end
A=diag(2*ones(1,n));
u=zeros(n-2,1);
lamda=zeros(n-1,1);
c=zeros(n,1);
fori=2:n-1
u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));
lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));
c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));
A(i,i+1)=u(i-1);
A(i,i-1)=lamda(i);
end
c(1)=2*s02;
c(n)=2*s10;
m=followup(A,c)
h=x(index+1)-x(index);
syms t;
f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+
y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/ h-m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h
f0=subs(f,'t',x0)
运行结果（方程S、弯矩M和插值函数f的值）为：
S0 = 0.9798652000
S1 = 0.9685577748
S2 = 0.9177710000
S3 = 0.8590474261
S4 = 0.8080348000
S5 = 0.7592534958
S6 = 0.6386092000
S7 = 0.4258081535
S8 =0.3843735000
f =125*((8825841099186633*t)/4503599627370496 - 8825841099186633/45035996273704960)*(t - 2/5)^2 - 125*((8266546267222895*t)/4503599627370496 - 8266546267222895/9007199254740992)*(t - 1/5)^2 + 25*((20271*t)/100000 - 20271/500000)*(t - 2/5)^2 - 25*((1321529499150801*t)/2251799813685248 - 1321529499150801/5629499534213120)*(t - 1/5)^2
3.画出y=S(x)的图形
（1）满足自然边界条件时
程序为：
x=0.2:0.2:1;
y=[0.9798652 0.9177710 0.8080348 0.6386093 0.3843735];
xi=0.2:0.01:1.0;
yi=interp1(x,y,xi,'variational');
plot(x,y,'o',xi,yi,'k-');
输出图形为：
图1 自然边界条件下的y=S(x)的图形
（2）满足第一类边界条件时
程序为：
x=[0.2 0.4 0.6 0.8 1.0 ];
y=[0.9798652 0.9177710 0.8080348 0.6386093 0.3843735];
pp=csape(x,y,'complete',[0.20271,1.55741]);
%complete代表一次插值
xi=0.2:0.01:1.0;
yi=ppval(pp,xi);
pp.coefsplot(x,y,'o',xi,yi);
输出图形为：
图2 第一类边界条件下的y=S(x)的图形
六、结果讨论和分析
在插值法中，拉格朗日插值法的公式结构整齐紧凑，在理论分析中十分方便，然而在计算中，当插值点增加或减少一个时，所对应的基本多项式就需要全部重新计算，于是整个公式都会变化，非常繁琐。此外，当插值点比较多的时候，拉格朗日插值多项式的次数可能会很高，因此具有数值不稳定的特点，也就是说尽管在已知的几个点取到给定的数值，但在附近却会和“实际上”的值之间有很大的偏差，这类现象也被称为龙格现象。分段线性插值，逼近程度虽然好，但光滑性差，即不能保证节点处插值函数的导数连续，从而不能满足某些工
程技术上的要求。分段三次Hermite插值，逼近程度好，光滑性也有所提高，但也增加了更多的条件，不太实用。而三次样条插值多项式，结合了二者的优点，即逼近程度好，光滑性强，不需要增加太多的条件，很实用。
实验过程中遇到了许多问题，首先就是MATLAB基础不太扎实，编程过程比较繁琐，程序比较复杂，在查阅大量资料并询问同学后，最终得以解决。其次，在编写满足第一类边界条件下的y=S(x)图形时，没有注意到interp1语句无法使用边界条件,致使编写出的程序出现错误，最后使用可以选择样条的边界条件的csape语句才得以解决问题。在工程上，构造三次样条插值函数通常有两种方法：一是以给定插值结点处得二阶导数值作为未知数来求解，而工程上称二阶导数为弯矩，因此，这种方法成为三弯矩插值。二是以给定插值结点处得一阶导数作为未知数来求解，而一阶导数又称为斜率，因此，这种方法称为三斜率插值。通过本次实验，我对三次样条插值的理解更加深刻，熟悉了通过求解三弯矩方程组得出曲线函数组的过程。同时也让我认识到三次样条插值的运用意义。
插值法数值上机实验报告
插值法数值上机实验报告 实验题目： 利用下列条件做插值逼近，并与R (x) 的图像比较 考虑函数：R x y=1 1+x2 （1）用等距节点X i=?5+i,i=0,1,...,10.给出它的10次Newton插值多项式的图像； π),i=0,1,...,20.给出它的20次Lagrange插值多项式（2）用节点X i=5cos(2i+1 42 的图像； （3）用等距节点X i=?5+i,i=0,1,...,10.给出它的分段线性插值函数的图像；（4）用等距节点X i=?5+i,i=0,1,...,10.给出它的三次自然样条插值函数的图像； （5）用等距节点X i=?5+i,i=0,1,...,10.给出它的分段三次Hermite插值函数的图像； 实验图像结果：
实验结果分析： 1.为了验证Range现象，我还特意做了10次牛顿插值多项式和20次牛顿插值多项式的对比图像，结果如下图（图对称，只截取一半） 可以看出，Range现象在高次时变得更加明显。这也是由于高次多项式在端点处的最值随次数的变大很明显。可以料定高次多项式在两侧端点处剧烈震荡，在更小的间距内急剧上升然后下降，Range现象非常明显。
2.分析实验（2）的结果，我们会惊讶地发现，由于取21个点逼近，原本预料的Range现象会很明显，但这里却和f(x)拟合的很好。（即下图中Lagrange p(x)的图像）。可是上图中取均匀节点的20次牛顿多项式逼近的效果在端点处却很差。料想是由于节点X i=5cos2i+1 42 π ,i=0,1,...,20 取得很好。由书上第五章的 知识，对于函数y=1 1+x ，y 1 2对应的cherbyshev多项式的根恰好为X i= 5cos2i+1 42 π ,i=0,1,...,20 。由于所学限制，未能深入分析。 （3）比较三次样条插值图像和Hermit插值图像对原函数图像的逼近情形。见下图：
样条插值函数与应用
样条插值函数及应用
摘要 样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。 在实际问题中所遇到许多函数往往很复杂，有些甚至是很难找到解析表达式的。根据函数已有的数据来计算函数在一些新的点处的函数值，就是插值法所需要解决的问题。 插值法是数值逼近的重要方法之一，它是根据给定的自变量值和函数值,求取未知函数的近似值。早在一千多年前,我国科学家就在研究历法时就用到了线性插值和二次插值。而在实际问题中,有许多插值函数的曲线要求具有较高的光滑性,在整个曲线中,曲线不但不能有拐点,而且曲率也不能有突变。因此，对于插值函数必须二次连续可微且不变号 ,这就需要用到三次样条插值。 关键词三次样条函数；插值法
目录 引言 0 第一章三次样条插值 (1) 1.1 样条插值函数简介 (1) 1.2 三次样条函数应用 (2) 第二章AMCM91A 估计水塔水流量 (4) 2.1 理论分析及计算 (5) 2.2运用MATLAB软件计算 (8) 参考文献 (13)
引言 样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。上世纪四十年代，在研究数据处理的问题中引出了样条函数，例如，在1946年Schoenberg将样条引入数学,即所谓的样条函数，直到五十年代，还多应用于统计数据的处理方面，从六十年代起，在航空、造船、汽车等行业中，开始大量采用样条函数。 在我国，从六十年代末开始，从船体数学放样到飞机外形设计，逐渐出现了一个使用样，逐渐出现了一个使用样条函数的热潮，并推广到数据处理的许多问题中。 在实际生活中有许多计算问题对插值函数的光滑性有较高的要求，例如飞机机翼外形、发动机进、排气口都要求有连续的二阶导数，用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时其函数值整体上能很好地逼近被插函数，相应的导数值也收敛于被插函数的导数值，不会发生“龙格现象”。 现在国内外学者对这方面的研究也越来越重视，根据我们的需要来解决不同的问题，而且函数的形式也在不断地改进，长期以来很多学者致力于样条插值的研究，对三次样条的研究已相当成熟。
三次样条插值、拉格朗日插值、herminte插值
三次样条插值： function s=spline(x0,y0,y2l,y2n,x) n=length(x0); km=length(x); a(1)=-0.5; b(1)=3*(y0(2)-y0(1))/(2*(x0(2)-x0(1))); for j=1:n-1 h(j)=x0(j+1)-x0(j); end for j=2:n-1 alpha(j)=h(j-1)/(h(j-1)+h(j)); beta(j)=3*((1-alpha(j))*y0(j)-y(j-1)/h(j-1)+alpha(j)*(y0(j+1)-y0(j))/h(j)); a(j)=-alpha(j)/(2+(1-alpha(j))*a(j-1)); b(j)=(beta(j)-(1-alpha(j))*b(j-1))/(2+(1-alpha(j))*a(j-1)); end m(n)=(3*(y0(n)-y0(n-1))/h(n-1)+y2n*h(n-1)/2-b(n-1))/(2+a(n-1)); for j=(n-1):-1:1 m(j)=a(j)*m(j+1)+b(j); end for k=1:km for j=1:(n-1) if ((x(k)>x0(j))&(x(k) 数值分析实验报告记录 ————————————————————————————————作者：————————————————————————————————日期： 数值分析实验报告 （第二章） 实验题目： 分别用二分法、牛顿迭代法、割线法、史蒂芬森迭代法求方程 的根，观察不同初始值下的收敛性，并给出结论。 问题分析： 题目有以下几点要求： 1.不同的迭代法计算根，并比较收敛性。 2.选定不同的初始值，比较收敛性。 实验原理： 各个迭代法简述 二分法：取有根区间的重点，确定新的有根区间的区间长度仅为区间长度的一版。对压缩了的有根区间重复以上过程，又得到新的有根区间，其区间长度为的一半，如此反复，……，可得一系列有根区间，区间收敛到一个点即为根。 牛顿迭代法：不动点迭代法的一种特例，具有局部二次收敛的特性。迭代格式为 割线法：是牛顿法的改进，具有超线性收敛的特性，收敛阶为1.618. 迭代格式为 史蒂芬森迭代法：采用不动点迭代进行预估校正。至少是平方收敛的。迭代格式为 这里可采用牛顿迭代法的迭代函数。 实验内容： 1.写出该问题的函数代码如下： function py= f(x) syms k; y=(k^2+1)*(k-1)^5; yy=diff(y,k); py(1)=subs(y,k,x); py(2)=subs(yy,k,x); end 2.分别写出各个迭代法的迭代函数代码如下： 二分法： function y=dichotomie(a,b,e) i=2; m(1)=a; while abs(a-b)>e t=(a+b)/2; s1=f(a); s2=f(b); s3=f(t); if s1(1)*s3(1)<=0 b=t; else a=t; end m(i)=t; i=i+1; end y=[t,i+1,m]; end 牛顿迭代法： function y=NewtonIterative(x,e) i=2; en=2*e;m(1)=x; while abs(en)>=e s=f(x); t=x-s(1)/s(2); en=t-x; x=t; m(i)=t; i=i+1; end y=[x,i+1,m]; end 牛顿割线法： function y=Secant(x1,x2,e) i=3; m(1)=x1,m(2)=x2; while abs(x2-x1)>=e s1=f(x1); s2=f(x2); t=x2-(x2-x1)*s2(1)/(s2(1)-s1( 1)); x1=x2; x2=t; m(i)=t; i=i+1; end 例1 设)(x f 为定义在[0,3]上的函数，有下列函数值表： 且2.0)('0=x f ，1)('3-=x f ，试求区间[0,3]上满足上述条件的三次样条插值函数)(x s 本算法求解出的三次样条插值函数将写成三弯矩方程的形式： ) ()6()() 6()(6)(6)(211123 13 1j j j j j j j j j j j j j j j j x x h h M y x x h h M y x x h M x x h M x s -- + -- + -+ -= +++++其中，方程中的系数 j j h M 6， j j h M 61+， j j j j h h M y )6(2- ， j j j 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) - 1 h(i) = IndVar(i + 1) - IndVar(i); end % 为向量μ赋值 mu = zeros(1, length(h)); for i = 1 : length(mu) - 1 mu(i) = h(i) / (h(i) + h(i + 1)); end mu(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)); 实验报告：牛顿差值多项式&三次样条 问题：在区间[-1,1]上分别取n=10、20用两组等距节点对龙格函数21()25f x x 作多项式插值及三次样条插值，对每个n 值，分别画出插值函数及()f x 的图形。 实验目的：通过编程实现牛顿插值方法和三次样条方法，加深对多项式插值的理解。应用所编程序解决实际算例。 实验要求： 1． 认真分析问题，深刻理解相关理论知识并能熟练应用； 2． 编写相关程序并进行实验； 3． 调试程序，得到最终结果； 4． 分析解释实验结果； 5． 按照要求完成实验报告。 实验原理： 详见《数值分析 第5版》第二章相关内容。 实验内容： （1）牛顿插值多项式 1.1 当n=10时： 在Matlab 下编写代码完成计算和画图。结果如下： 代码： clear all clc x1=-1:0.2:1; y1=1./(1+25.*x1.^2); n=length(x1); f=y1(:); for j=2:n for i=n:-1:j f(i)=(f(i)-f(i-1))/(x1(i)-x1(i-j+1)); end end syms F x p ; F(1)=1;p(1)=y1(1); for i=2:n F(i)=F(i-1)*(x-x1(i-1)); p(i)=f(i)*F(i); end syms P P=sum(p); P10=vpa(expand(P),5); x0=-1:0.001:1; y0=subs(P,x,x0); y2=subs(1/(1+25*x^2),x,x0); plot(x0,y0,x0,y2) grid on xlabel('x') ylabel('y') P10即我们所求的牛顿插值多项式，其结果为：P10（x）=-220.94*x^10+494.91*x^8-9.5065e-14*x^7-381.43*x^6-8.504e-14*x^5+123.36*x^4+2.0 202e-14*x^3-16.855*x^2-6.6594e-16*x+1.0 并且这里也能得到该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.1）。 Fig.1 牛顿插值多项式（n=10）函数和原函数图形 从图形中我们可以明显的观察出插值函数在两端点处发生了剧烈的波动，产生了极大的误差，即龙格现象，当n=20时，这一现象将更加明显。 1.2 当n=20时： 对n=10的代码进行修改就可以得到n=20时的代码。将“x1=-1:0.2:1;”改为“x1=-1:0.1:1;”即可。运行程序，我们得到n=20时的牛顿插值多项式，结果为：P20(x)= 260188.0*x^20 - 1.0121e6*x^18 + 2.6193e-12*x^17 + 1.6392e6*x^16 + 2.248e-11*x^15 - 1.4429e6*x^14 - 4.6331e-11*x^13 + 757299.0*x^12 + 1.7687e-11*x^11 - 245255.0*x^10 + 2.1019e-11*x^9 + 49318.0*x^8 + 3.5903e-12*x^7 - 6119.2*x^6 - 1.5935e-12*x^5 + 470.85*x^4 + 1.3597e-14*x^3 - 24.143*x^2 - 1.738e-14*x + 1.0 同样的，这里得到了该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.2）。 . 牛顿插值法一、实验目的：学会牛顿插值法，并应用算法于实际问题。 x?x)f(二、实验内容：给定函数，已知： 4832401.2)?.?1449138f(2.f.f(20)?1.414214(2.1) 549193.)?1f(2.4516575(f2.3)?1. 三、实验要求：以此作为函数2.15插值多项式在处的值，用牛顿插值法求4 次Newton（ 1）2.15?N(2.15)。在MATLAB中用内部函数ezplot绘制出的近似值4次Newton插值多项式的函数图形。 （2）在MATLAB中用内部函数ezplot可直接绘制出以上函数的图形，并与作出的4次Newton插值多项式的图形进行比较。 四、实验过程： 1、编写主函数。打开Editor编辑器，输入Newton插值法主程序语句： function [y,L]=newdscg(X,Y,x) n=length(X); z=x; A=zeros(n,n);A(:,1)=Y';s=0.0; p=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end end C=A(n,n); for k=(n-1):-1:1 C=conv(C,poly(X(k))); d=length(C);C(d)=C(d)+A(k,k); end y(k)= polyval(C, z); L(k,:)=poly2sym(C); 0 / 3 . %%%%%%%%%%%%%%%%%% t=[2,2.1,2.2,2.3,2.4]; fx=sqrt(t); wucha=fx-Y; 以文件名newdscg.m保存。 2、运行程序。 （1）在MATLAB命令窗口输入： >> X=[2,2.1,2.2,2.3,2.4]; Y =[1.414214,1.449138,1.483240,1.516575,1.549193]; x=2.15;[y,P]=newdscg(X,Y,x) 回车得到： natural cubic spine二阶导数求解 1 Natural cubic spine is The formula (3.37) in book ‘Numerical Recipes in Fortran 77’ is: (1) Where . Formula (1) can be simplified as: (2) Where: (3) Formula (2) can be written as: (4) This is a ‘m=n-2’ dimensional equation, then change the subscript. Take: , and (4) can be changed into: (5) To solve equation (5), we can get , then we can get . 三次样条插值法使用时，一般要对y均分取值。 2 The second derivative is known: , (6) When : (7) When : (8) So formula (4) can be written as: (9) Then we can solve equation (9). 3 The first derivative is known. The formula (3.35) is : (6) Where: Case 1: , (6) can be changed into: (7) And , (7) can be written as: (8) (9) Case 2: , (6) can be written as: (10) And , (10) can be written as: (11) (12) , (13) j=2: 实验二插值法 1、实验目的： 1、掌握直接利用拉格郎日插值多项式计算函数在已知点的函数值；观察拉格郎日插值的龙格现象。 2、了解Hermite插值法、三次样条插值法原理，结合计算公式，确定函数值。 2、实验要求: 1)认真分析题目的条件和要求，复习相关的理论知识，选择适当的解决方案和算法； 2)编写上机实验程序，作好上机前的准备工作； 3)上机调试程序，并试算各种方案，记录计算的结果（包括必要的中间结果）； 4)分析和解释计算结果； 5)按照要求书写实验报告； 3、实验内容： 1) 用拉格郎日插值公式确定函数值；对函数f(x)进行拉格郎日插值，并对f(x)与插值多项式的曲线作比较。 已知函数表：（0.56160,0.82741）、（0.56280,0.82659）、（0.56401,0.82577）、（0.56521,0.82495）用三次拉格朗日插值多项式求x=0.5635时函数近似值。 2) 求满足插值条件的插值多项式及余项 1) 4、题目：插值法 5、原理： 拉格郎日插值原理： n次拉格朗日插值多项式为：L n (x)=y l (x)+y 1 l 1 (x)+y 2 l 2 (x)+…+y n l n (x) n=1时，称为线性插值， L 1(x)=y (x-x 1 )/(x -x 1 )+y 1 (x-x )/(x 1 -x )=y +(y 1 -x )(x-x )/(x 1 -x ) n=2时，称为二次插值或抛物线插值， L 2(x)=y (x-x 1 )(x-x 2 )/(x -x 1 )/(x -x 2 )+y 1 (x-x )(x-x 2 )/(x 1 -x )/(x 1 -x 2 )+y 2 (x -x 0)(x-x 1 )/(x 2 -x )/(x 2 -x 1 ) n=i时， Li= （X-X0）……(X-X i-1)(x-x i+1) ……(x-x n) （X-X0）……(X-X i-1)(x-x i+1) ……(x-x n) 6、设计思想： 拉格朗日插值法是根据n + 1个点x0, x1, ... x n(x0 < x1 < ... x n)的函数值f (x0), f (x1) , ... , f (x n)推出n次多項式p(x)，然后n次多項式p (x)求出任意的点x对应的函数值f (x)的算法。 7、对应程序： 1 ) 三次拉格朗日插值多项式求x=0.5635时函数近似值 #include"stdio.h" #define n 5 void main() { int i,j; float x[n],y[n]; float x1; float a=1; float b=1; float lx=0; printf("\n请输入想要求解的X：\n x="); scanf("%f",&x1); printf("请输入所有点的横纵坐标：\n"); for(i=1;i 实验报告：牛顿差值多项式&三次样条 问题：在区间[-1,1]上分别取n=10、20用两组等距节点对龙格函数2 1()25f x x 作多项式插值及三次样条插值，对每个n 值，分别画出插值函数及()f x 的图形。 实验目的：通过编程实现牛顿插值方法和三次样条方法，加深对多项式插值的理解。应用所编程序解决实际算例。 实验要求： 1． 认真分析问题，深刻理解相关理论知识并能熟练应用； 2． 编写相关程序并进行实验； 3． 调试程序，得到最终结果； 4． 分析解释实验结果； 5． 按照要求完成实验报告。 实验原理： 详见《数值分析 第5版》第二章相关内容。 实验内容： （1）牛顿插值多项式 1.1 当n=10时： 在Matlab 下编写代码完成计算和画图。结果如下： 代码： clear all clc x1=-1:0.2:1; y1=1./(1+25.*x1.^2); n=length(x1); f=y1(:); for j=2:n for i=n:-1:j f(i)=(f(i)-f(i-1))/(x1(i)-x1(i-j+1)); end end syms F x p ; F(1)=1;p(1)=y1(1); for i=2:n F(i)=F(i-1)*(x-x1(i-1)); p(i)=f(i)*F(i); end syms P P=sum(p); P10=vpa(expand(P),5); x0=-1:0.001:1; y0=subs(P,x,x0); y2=subs(1/(1+25*x^2),x,x0); plot(x0,y0,x0,y2) grid on xlabel('x') ylabel('y') P10即我们所求的牛顿插值多项式，其结果为：P10（x）=-220.94*x^10+494.91*x^8-9.5065e-14*x^7-381.43*x^6-8.504e-14*x^5+123.36*x^4+2.0202e-1 4*x^3-16.855*x^2-6.6594e-16*x+1.0 并且这里也能得到该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.1）。 Fig.1 牛顿插值多项式（n=10）函数和原函数图形 从图形中我们可以明显的观察出插值函数在两端点处发生了剧烈的波动，产生了极大的误差，即龙格现象，当n=20时，这一现象将更加明显。 1.2 当n=20时： 对n=10的代码进行修改就可以得到n=20时的代码。将“x1=-1:0.2:1;”改为“x1=-1:0.1:1;”即可。运行程序，我们得到n=20时的牛顿插值多项式，结果为：P20(x)= 260188.0*x^20 - 1.0121e6*x^18 + 2.6193e-12*x^17 + 1.6392e6*x^16 + 2.248e-11*x^15 - 1.4429e6*x^14 - 4.6331e-11*x^13 + 757299.0*x^12 + 1.7687e-11*x^11 - 245255.0*x^10 + 2.1019e-11*x^9 + 49318.0*x^8 + 3.5903e-12*x^7 - 6119.2*x^6 - 1.5935e-12*x^5 + 470.85*x^4 + 1.3597e-14*x^3 - 24.143*x^2 - 1.738e-14*x + 1.0 同样的，这里得到了该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.2）。 CENTRAL SOUTH UNIVERSITY 数值分析实验报告 三次样条插值方法的应用 一、问题背景 分段低次插值函数往往具有很好的收敛性，计算过程简单，稳定性好，并且易于在在电子计算机上实现，但其光滑性较差，对于像高速飞机的机翼形线船体放样等型值线往往要求具有二阶光滑度，即有二阶连续导数，早期工程师制图时，把富有弹性的细长木条（即所谓的样条）用压铁固定在样点上，在其他地方让他自由弯曲，然后沿木条画下曲线，称为样条曲线。样条曲线实际上是由分段三次曲线并接而成，在连接点即样点上要求二阶导数连续，从数学上加以概括就得到数学样条这一概念。下面我们讨论最常用的三次样条函数及其应用。 二、数学模型 样条函数可以给出光滑的插值曲线（面），因此在数值逼近、常微分方程和偏微分方程的数值解及科学和工程的计算中起着重要的作用。 设区间[]b ,a 上给定有关划分b x x n =<<<= 10x a ，S 为[]b ,a 上满足下面条件的函数。 ● )(b a C S ,2∈； ● S 在每个子区间[]1,+i i x x 上是三次多项式。 则称S 为关于划分的三次样条函数。常用的三次样条函数的边界条件有三种类型： ● Ⅰ型 ()()n n n f x S f x S ''0'',==。 ● Ⅱ型 ()()n n n f x S f x S ''''0'''',==，其特殊情况为()()0''''==n n x S x S 。 ● Ⅲ型 ()() 3,2,1,0,0==j x S x S n j j ，此条件称为周期样条函数。 鉴于Ⅱ型三次样条插值函数在实际应用中的重要地位，在此主要对它进行详细介绍。 三、算法及流程 按照传统的编程方法，可将公式直接转换为MATLAB 可是别的语言即可；另一种是运用矩阵运算，发挥MATLAB 在矩阵运算上的优势。两种方法都可以方便地得到结果。方法二更直观，但计算系数时要特别注意。这里计算的是方法一的程序，采用的是Ⅱ型边界条件，取名为spline2.m 。 Matlab 代码如下： function s=spline2(x0,y0,y21,y2n,x) %s=spline2(x0,y0,y21,y2n,x) %x0,y0 are existed points,x are insert points,y21,y2n are the second 数值计算方法作业 实验4.3 三次样条差值函数 实验目的： 掌握三次样条插值函数的三弯矩方法。 实验函数: dt e x f x t ? ∞ -- = 2 221)(π 实验内容： （1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值； （3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线 比较插值结果。 实验4.5 三次样条差值函数的收敛性 实验目的： 多项式插值不一定是收敛的，即插值的节点多，效果不一定好。对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。 实验内容： 按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。 实验要求： （1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情 况，分析所得结果并与拉格朗日插值多项式比较； （2） 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子，考 虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一 算法描述： 拉格朗日插值： 错误!未找到引用源。 其中错误!未找到引用源。是拉格朗日基函数，其表达式为：() ∏ ≠=--=n i j j j i j i 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 j i j i j i 三样条插值： 所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a 实验报告内容： 一：不动点迭代法解方程 二：牛顿插值法的MATLAB实现 完成日期：2012年6月21日星期四 数学实验报告一 日期：2012-6-21 所以，确定初值为x0=1 二：不断迭代 算法： 第一步：将f(x0)赋值给x1 第二步：确定x1-x0的绝对值大小，若小于给定的误差值，则将x1当做方程的解，否则回到第一步 编写计算机程序： clear f=inline('0.5*sin(x)+0.4'); x0=1; x1=f(x0); k=1; while abs(x1-x0)>=1.0e-6 x0=x1; x1=f(x0); k=k+1; fprintf('k=%.0f,x0=%.9f,x1=%.9f\n',k,x0,x1) end 显示结果如下： k=2,x0=0.820735492,x1=0.765823700 k=3,x0=0.765823700,x1=0.746565483 k=4,x0=0.746565483,x1=0.739560873 k=6,x0=0.736981783,x1=0.736027993 k=7,x0=0.736027993,x1=0.735674699 k=8,x0=0.735674699,x1=0.735543758 k=9,x0=0.735543758,x1=0.735495216 k=10,x0=0.735495216,x1=0.735477220 k=11,x0=0.735477220,x1=0.735470548 k=12,x0=0.735470548,x1=0.735468074 k=13,x0=0.735468074,x1=0.735467157 >>。。。 以下是程序运行截图： 学习报告—— 三次样条函数插值问题的讨论 班级：数学二班 学号：152111033 姓名：刘楠楠 样条函数： 由一些按照某种光滑条件分段拼接起来的多项式组成的函数；最常用的样条函数为三次样条函数，即由三次多项式组成，满足处处有二阶连续导数。 一、三次样条函数的定义： 对插值区间[,]a b 进行划分，设节点011n n a x x x x b -=<< <<=，若 函数2()[,]s x c a b ∈在每个小区间1[,]i i x x +上是三次多项式，则称其为三次样条函数。如果同时满足()()i i s x f x = (0,1,2)i n =，则称()s x 为()f x 在 [,]a b 上的三次样条函数。 二、三次样条函数的确定： 由定义可设：101212 1(),[,] (),[,]()(),[,] n n n s x x x x s x x x x s x s x x x x -∈??∈?=???∈?其中()k s x 为1[,]k k x x -上的三次 多项式，且满足11(),()k k k k k k s x y s x y --== (1,2,,k n = 由2()[,]s x C a b ∈可得：''''''()(),()(),k k k k s x s x s x s x -+-+== 有''1()(),k k k k s x s x -++= ''''1()(),(1 ,2,,1)k k k k s x s x k n -+ +==-， 已知每个()k s x 均为三次多项式，有四个待定系数，所以共有4n 个待定系数，需要4n 个方程才能求解。前面已经得到22(1)42n n n +-=-个方程，因此要唯一确定三次插值函数，还要附加2个条件，一般上，实际问题通常对样条函数在端点处的状态有要求，即所谓的边界条件。 1、第一类边界条件：给定函数在端点处的一阶导数，即 ''''00(),()n n s x f s x f == 2、第二类边界条件：给定函数在端点处的二阶导数，即 《数值分析》课程实验报告 用拉格朗日和牛顿插值法求解函数值 算法名称用拉格朗日和牛顿插值法求函数值 学科专业xxxxx 作者姓名xxxx 作者学号xxxxx 作者班级xxxxxx xxx大学 二〇一五年十二月 《数值分析》课程实验报告 得到的近似值为。 拉格朗日插值模型简单，结构紧凑，是经典的插值法。但是由于拉格朗日的插值多项式和每个节点都有关，当改变节点个数时，需要重新计算。且当增大插值阶数时容易出现龙格现象。 2.牛顿插值法 在命令窗口输入： x=[ ]; y=[ ]; xt=; [yt,N]=NewtInterp(x,y,xt) z=::2; yz=subs(N,'t',z); figure; plot(z,sqrt(z),'--r',z,yz,'-b') hold on plot(x,y,'marker','+') hold on plot(xt,yt,'marker','o') h=legend('$\sqrt{x}$','牛顿','$(x_k,y_k)$','$x=$'); set(h,'Interpreter','latex') xlabel('x') ylabel('y') 得到结果及图像如下： yt = N = - *t^4 + *t^3 - *t^2 + *t + 得到√的近似值为，插值函数为 N =- *t^4 + *t^3 - *t^2 + *t + ， 其计算精度是相当高的。 Lagrange插值法和Newton插值法解决实际问题中关于只提供复杂的离散数据的函数求值问题，通过将所考察的函数简单化，构造关于离散数据实际函数f(x)的近似函数P(x)，从而可以计算未知点出的函数值，是插值法的基本思路。 实际上Lagrange插值法和Newton插值法是同一种方法的两种变形，其构造拟合函数的思路是相同的，而实验中两个实际问题用两种算法计算出结果是相同的。 牛顿插值法 一、实验目的：学会牛顿插值法，并应用算法于实际问题。 二、实验内容：给定函数 x x f =)(，已知： 414214.1)0.2(=f 449138.1)1.2(=f 483240.1)2.2(=f 516575.1)3.2(=f 549193.1)4.2(=f 三、实验要求： （1）用牛顿插值法求4次Newton 插值多项式在2.15处的值，以此作为函数的近似值)15.2(15.2N ≈。在MATLAB 中用内部函数ezplot 绘制出4次Newton 插值多项式的函数图形。 （2）在MATLAB 中用内部函数ezplot 可直接绘制出以上函数的图形，并与作出的4次Newton 插值多项式的图形进行比较。 四、实验过程： 1、编写主函数。打开Editor 编辑器，输入Newton 插值法主程序语句： function [y,L]=newdscg(X,Y,x) n=length(X); z=x; A=zeros(n,n);A(:,1)=Y';s=0.0; p=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end end C=A(n,n); for k=(n-1):-1:1 C=conv(C,poly(X(k))); d=length(C);C(d)=C(d)+A(k,k); end y(k)= polyval(C, z); L(k,:)=poly2sym(C); %%%%%%%%%%%%%%%%%% t=[2,2.1,2.2,2.3,2.4]; fx=sqrt(t); wucha=fx-Y; 以文件名newdscg.m保存。 2、运行程序。 （1）在MATLAB命令窗口输入： >> X=[2,2.1,2.2,2.3,2.4]; Y =[1.414214,1.449138,1.483240,1.516575,1.549193]; x=2.15;[y,P]=newdscg(X,Y,x) 回车得到： y =1.4663 wucha =1.0e-06 * -0.4376 -0.3254 -0.3026 0.0888 0.3385 P = - (4803839603609061*x^4)/2305843009213693952 + (7806239355294329*x^3)/288230376151711744 - (176292469178709*x^2)/1125899906842624 + (1624739243112817*x)/2251799813685248 + 1865116246031207/4503599627370496 （2）在MATLAB命令窗口输入： >> v=[0,6,-1,3]; >> ezplot(P),axis(v),grid >> hold on >> x=0:0.1:6; >> yt=sqrt(x);plot(x,yt,':') >> legend('插值效果','原函数') >> xlabel('X') >> ylabel('Y') >>title('Newton插值与原函数比较') 回车即可得到图像1-1。 计算方法实验 题目： 班级： 学号： 姓名： 目录 计算方法实验 (1) 1 实验目的 (3) 2 实验步骤 (3) 2.1环境配置： (3) 2.2添加头文件 (3) 2.3主要模块 (3) 3 代码 (4) 3.1主程序部分 (4) 3.2多项式方程部分 (4) 3.3核心算法部分 (8) 3.4数据结构部分 (13) 4运行结果 (19) 4.1拉格朗日插值法运行结果 (19) 4.2牛顿插值法运行结果 (20) 4.3多项式拟合运行结果 (20) 5总结 (21) 拉格朗日插值法 (21) 牛顿插值法 (21) 多项式拟合 (21) 6参考资料 (22) 1 实验目的 1.通过编程对拉格朗日插值法、牛顿插值法以及多项式拟合数据的理解 2.观察上述方法的计算稳定性和求解精度并比较各种方法利弊 2 实验步骤 2.1环境配置： VS2013，C++控制台程序 2.2添加头文件 #include "stdio.h" #include "stdlib.h" #include "stdafx.h" 2.3主要模块 程序一共分成三层，最底层是数据结构部分，负责存储数据，第二层是交互部分，即多项式方程部分，负责输入输出获得数据，最上层是核心的算法部分，负责处理已获得的数据。具体功能如下： ●数据结构部分 数据结构部分是整个程序的最底层，负责存储部分。因方程系数作为数据元素插入和删除操作较少，而顺序表空间利用率大且查看方便，故此程序选用顺序表保存系数。数据结构文件中写的是有关顺序表的所有基本操作以供其他文件调用。本次实验使用列主元高斯消元法作为求解方程组的方法，所以也用了二维顺序表存储数组。综上，数据结构部分文件是前两个试验的文件内容和，稍作修改。 ●常系数微分方程部分 多项式方程部分是程序的第二层，内容主要是常系数微分方程导数的计算和显示菜单部分。 ●算法部分 算法部分分为两个文件，一个是插值部分，一个是拟合部分。 插值部分文件负责有关插值的核心算法，处于整个程序最上层部分，负责拉格朗日插值法和牛顿插值法的具体实现过程。调用方程文件的函数，将获得的数据进行处理运算，将结果返回给方程主函数和输出的第二层。每种方法有两个函数，一个为仅仅实现一次插值的算法，另一个是和方程部分联系的数值分析实验报告记录
三次样条插值课后题集
数值分析实验报告-插值、三次样条Word版
牛顿插值法试验报告
三次样条插值法
插值法实验报告
数值分析实验报告-插值、三次样条(教育教学)
三次样条插值方法的应用
数值分析作业-三次样条插值
matlab(迭代法-牛顿插值)Word版
关于三次样条插值函数的学习报告(研究生)资料
数值分析课程实验报告-拉格朗日和牛顿插值法
牛顿插值法实验报告
计算方法上机作业插值与拟合实验报告