数值计算方法作业
实验 三次样条差值函数
实验目的:
掌握三次样条插值函数的三弯矩方法。
实验函数:
dt e
x f x
t ?
∞
--
=
2
221)(π
实验内容:
(1) 编程实现求三次样条插值函数的算法,分别考虑不同的边界条件; (2) 计算各插值节点的弯矩值;
(3) 在同一坐标系中绘制函数f(x),插值多项式,三次样条插值多项式的曲线
比较插值结果。
实验 三次样条差值函数的收敛性
实验目的:
多项式插值不一定是收敛的,即插值的节点多,效果不一定好。对三次样条插值函数如何呢理论上证明三次样条插值函数的收敛性是比较困难的,通过本实验可以证明这一理论结果。
实验内容:
按照一定的规则分别选择等距或非等距的插值节点,并不断增加插值节点的个数。
实验要求:
(1) 随着节点个数的增加,比较被逼近函数和三样条插值函数的误差变化情
况,分析所得结果并与拉格朗日插值多项式比较;
(2) 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子,考
虑如下例子:某汽车制造商根据三次样条插值函数设计车门曲线,其中一
k
x
012345678910 k
y
k
y'
算法描述:
拉格朗日插值:
其中是拉格朗日基函数,其表达式为:
()
∏
≠
=
-
-
=
n
i
j
j j
i
j
i x
x
x
x
x
l
)
(
)
(
牛顿插值:
)
)...(
)(
](
,...
,
,
[
....
)
)(
](
,
,
[
)0
](
,
[
)
(
)
(
1
1
2
1
1
2
1
1
-
-
-
-
+
+
-
-
+
-
+
=
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
其中
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
-
-
=
-
-
=
-
-
=
-
)
/(
])
,...
,
[
]
,...
,
[
(
]
...
,
[
.
.
]
,
[
]
,
[
]
,
,
[
)
(
)
(
]
,
[
1
1
2
1
1
x
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 ] , [ ), 6 ( ) 6 ( ] 6 ) ( [ 6 ) ( 6 ) ( ) ( 1 1 1 1 1 3 1 3 1 i i i i i i i i i i i i i i i i i i i i i x x x h y M 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)(61111111 11 i i i i i i i i i i i i i i i i i i i i x 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 ='='='=')(,)(00 0 (2) 给定区间两端点的二阶导数M0,Mn,即n n n M y x S M y x S =''=''=''='')(,)(00 0 (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 221 100μλ 其中 n n n n n n n M u x x f m h d M m x x f h d )1(2]),[(6)1(2)],[(610000101 0-+-= -+-= -μλλ 那么解就可以为 ????????? ???????????=????????????? ????????????????? ??----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 n n n n x x f x x f h h d h h h h h h --+=+=+= μλ,那么解就可以为: ??????????? ?? ???????=?????????????????????????????? ??-------1221012101221100...2.............2..............................2..,,.......,..22 n n n n n n n d d d d d M M M M M n μλλμλμμλ 程序代码: 1拉格朗日插值函数 function f=lang(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标 %xi 插值点处的横坐标 %f 求得的拉格朗日插值多项式的值 n=length(X); f=0; for i=1:n l=1; for j=1:i-1 l=l.*(xi-X(j))/(X(i)-X(j)); end ; for j=i+1:n l=l.*(xi-X(j))/(X(i)-X(j)); end ;%拉格朗日基函数 f=f+l*Y(i); end fprintf('%d\n',f) return 2 牛顿插值函数 function f=newton(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标 %xi 插值点处的横坐标 %f 求得的拉格朗日插值多项式的值 n=length(X); newt=[X',Y']; %计算差商表 for j=2:n for i=n:-1:1 if i>=j Y(i)=(Y(i)-Y(i-1))/(X(i)-X(i-j+1)); else Y(i)=0; end end newt=[newt,Y']; end %计算牛顿插值 f=newt(1,2); for i=2:n z=1; for k=1:i-1 z=(xi-X(k))*z; end f=f+newt(i-1,i)*z; end fprintf('%d\n',f) return 3三次样条插值第一类边界条件 function 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); end for i=2:n%求函数的二阶差商 f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end d(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-1 B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=1; A(n+1,n)=1; for i=1:n+1 A(i,i)=2; end for i=2:n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\d; syms x; for i=1:n Sx(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));%三样条插值函数表达式 end for i=1:n disp('S(x)='); fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1)); end for i=1:n if 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; end end disp('xi S'); fprintf('%d,%d\n',xi,S); return 4 三次样条插值第二类边界条件 function [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); end for i=2:n%求二阶差商 f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end d(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-1 B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=0; A(n+1,n)=0; for i=1:n+1 A(i,i)=2; end for i=2:n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\d; syms x; for i=1:n Sx(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)); end for i=1:n disp('S(x)='); fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1)); end for i=1:n if 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; end end disp('xi S'); fprintf('%d,%d\n',xi,S); return 5插值节点处的插值结果 clear clc X=[,,,,]; Y=[,,,,]; xi=; %xi=; disp('xi='); %disp('xi='); disp('拉格朗日插值结果'); lang(X,Y,xi); disp('牛顿插值结果'); newton(X,Y,xi); disp('三次样条第一类边界条件插值结果'); Threch1(X,Y,,,xi);%,分别为两端点处的一阶导数 disp('三次样条第二类边界条件插值结果'); Threch2(X,Y,0,,xi);%0,分别为两端点处的二阶导数 6将多种插值函数即原函数图像画在同一张图上 clear clc X=[,,,,]; Y=[,,,,]; a=linspace(0,,21); NUM=21; L=zeros(1,NUM); N=zeros(1,NUM); S=zeros(1,NUM); B=zeros(1,NUM); for i=1:NUM xi=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,,,xi);%三次样条函数第一类边界条件end plot(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 off 7增加插值节点观察误差变化 clear; clc; N=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:101 L(j)=lang(t,ft,val(j)); S(j)=Threch1(t,ft,,,val(j));%三样条第一类边界条件插值 end plot(a,Ini,'k')%原函数图象 hold on plot(val,L,'r')%拉格朗日插值函数图像 hold on plot(val,S,'b')%三次样条插值函数图像 str=sprintf('插值节点为%d时的插值效果',N); title(str); legend('原函数','拉格朗日插值','三次样条插值');%显示图例hold off figure end 8车门曲线 clear clc X=[0,1,2,3,4,5,6,7,8,9,10]; Y=[,,,,,,,,,,]; 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:nh(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i))/h(i); end for i=2:nf2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end d(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-1 B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=1; A(n+1,n)=1; for i=1:n+1 A(i,i)=2; end for i=2:n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\d; x=zeros(1,n); S=zeros(1,n); for i=1:n x(i)=X(i)+; 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; end plot(X,Y,'k'); hold on; plot(x,S,'o'); title('三次样条插值效果图'); legend('已知插值节点','三次样条插值'); hold off 实验结果: 1计算插值节点处的函数值 xi=时 Xi=时 2将多种插值函数即原函数图像画在同一张图上 增加插值节点观察误差变化 从上面三张图可以看出增加插值节点并不能改善差之效果车门曲线 样条插值函数及应用 摘要 样条函数具有广泛的应用,是现代函数论的一个十分活跃的分支,是计算方法的主要基础和工具之一,由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要,人们便更多地应用这个工具,也更深刻的认识了它的本质。 在实际问题中所遇到许多函数往往很复杂,有些甚至是很难找到解析表达式的。根据函数已有的数据来计算函数在一些新的点处的函数值,就是插值法所需要解决的问题。 插值法是数值逼近的重要方法之一,它是根据给定的自变量值和函数值,求取未知函数的近似值。早在一千多年前,我国科学家就在研究历法时就用到了线性插值和二次插值。而在实际问题中,有许多插值函数的曲线要求具有较高的光滑性,在整个曲线中,曲线不但不能有拐点,而且曲率也不能有突变。因此,对于插值函数必须二次连续可微且不变号 ,这就需要用到三次样条插值。 关键词三次样条函数;插值法 目录 引言 0 第一章三次样条插值 (1) 1.1 样条插值函数简介 (1) 1.2 三次样条函数应用 (2) 第二章AMCM91A 估计水塔水流量 (4) 2.1 理论分析及计算 (5) 2.2运用MATLAB软件计算 (8) 参考文献 (13) 引言 样条函数具有广泛的应用,是现代函数论的一个十分活跃的分支,是计算方法的主要基础和工具之一,由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要,人们便更多地应用这个工具,也更深刻的认识了它的本质。上世纪四十年代,在研究数据处理的问题中引出了样条函数,例如,在1946年Schoenberg将样条引入数学,即所谓的样条函数,直到五十年代,还多应用于统计数据的处理方面,从六十年代起,在航空、造船、汽车等行业中,开始大量采用样条函数。 在我国,从六十年代末开始,从船体数学放样到飞机外形设计,逐渐出现了一个使用样,逐渐出现了一个使用样条函数的热潮,并推广到数据处理的许多问题中。 在实际生活中有许多计算问题对插值函数的光滑性有较高的要求,例如飞机机翼外形、发动机进、排气口都要求有连续的二阶导数,用三次样条绘制的曲线不仅有很好的光滑度,而且当节点逐渐加密时其函数值整体上能很好地逼近被插函数,相应的导数值也收敛于被插函数的导数值,不会发生“龙格现象”。 现在国内外学者对这方面的研究也越来越重视,根据我们的需要来解决不同的问题,而且函数的形式也在不断地改进,长期以来很多学者致力于样条插值的研究,对三次样条的研究已相当成熟。 对样条函数及其插值问题的一点认识 样条函数是计算数学以及计算机辅助设计几何设计的重要工具。1946年,I. J. Schoenberg 著名的关于一元样条函数的奠定性论文“Contribution to the problem of application of equidistant data by analytic functions ”发表,建立了一元样条函数的理论基础。自此以后,关于样条函数的研究工作逐渐深入。随着电子计算机技术的不断进步,样条函数的理论以及应用研究得到迅速的发展和广泛的应用。经过数学工作者的努力,已经形成了较为系统的理论体系。 所谓(多项式)样条函数,乃指具有一定光滑性的分段(分片)多项式。一元n 次且n -1阶连续可微的样条函数具有如下的表示式: 1()()()()N n n j j j s x p x c x x x +==+--∞<<+∞∑[] 011,00,01,,...,,(1),...,(),,...,,n n n n N n N N u un u u u u x x x x x S x x x x ++++ +≥??=???-- 其中()n p x 是一元n 次多项式,12,,...,n x x x 称为样条节点。 用 n u +记u 的n 次截断幂: ,00,0u u u u +≥??=??? 则1,,...,,(1),...,()n n n N x x x x x ++--构成了由这样的样条函数全体组成的线形空间 []011,,...,,n N N S x x x x +的一组基。 1953年,I. J. Schoenberg 与A. Whitney 获得了判断一元样条插值结点组是否为适定结点组的准则,即著名的Schoenberg-Whitney 定理。1966年,H. B. Curry 与I. J. Schoenberg 引入一元B-样条,给出了重要的样条函数的B-样条基的表示方法。有了完整的B-样条基理论之后,样条函数逼近无论在理论研究还是在应用问题探讨方面都更加方便。此后,关于样条函数的理论以及应用的研究不断取得进展。特别地,随着计算机技术的飞速发展,人们进一步认识到样条如同多项式一般的计算方便性以及强于多项式的局部可调的灵活性、易存储性等诸多性质的重要意义,并把它应用到与科学计算相关的许多领域,比如数值逼近、微分方程数值解、计算机辅助几何设计、小波及有限元等。同时,样条函数在各个方面的推广也成为数学工作者们密切关注并开展积极研究的重要课题。 样条在多元方面的推广自1960年Birkhoff 与Garabedian 开始,但是,由于它的复杂性,这方面的研究工作不如一元样条函数那样顺利。1962年,C. deBoor 研究并证明了一些双三次内插样条的存在与唯一性。但其方法本质上只是一元样条函数的简单推广。1975年,王仁宏教授采用函数论与代数几何的理论,提出并运用光滑余因子协调法,建立了任意分割部分下的多元样条函数的基本理论框架。他引入了分片代数曲线概念,证明了多元样条的一个插值结点组为适定结点组的冲要条件是这些结点不同时位于一条非零分片代数曲线上。采用 《数值分析课程设计-三次样条插值》 报告 掌握三次样条插值函数的构造方法,体会三次样条插值函数对被逼近函数的近似。 三次样条插值函数边界条件由实际问题对三次样条插值在端点的状态要求给出。 以第1 边界条件为例,用节点处二阶导数表示三次样条插值函数,用追赶法求解相关方程组。通过Matlab 编制三次样条函数的通用程序,可直接显示各区间段三次样条函数体表达式,计算出已给点插值并显示各区间分段曲线图。 引言 分段低次样条插值虽然计算简单、稳定性好、收敛性有保证且易在电子计算机上实现,但只能保证各小段曲线在连接处的连续性,不能保证整件曲线的光滑性。利用样条插值,既可保持分段低次插值多项式,又可提高插值函数光滑性。故给出分段三次样条插值的构造过程算法步骤,利用Matlab软件编写三次样条插值函数通用程序,并通过数值算例证明程序的正确性。 三次样条函数的定义及特征 定义:设[a,b] 上有插值节点,a=x1<x2<…xn=b,对应函数值为y1,y2,?yn。若函数S(x) 满 足S(xj) = yj ( j = 1,2, ?,n ), S(x) 在[xj,xj+1] ( j =1,2,?,n-1)上都是不高于三的多项式(为了与其对应j 从1 开始,在Matlab 中元素脚标从1 开始)。当S(x) 在 [a,b] 具有二阶连续导 数。则称S(x) 为三次样条插值函数。要求S(x) 只需在每个子区 间[xj,xj+1] 上确定 1 个三次多项式,设为: Sj(x)=ajx3+bjx2+cjx+dj, (j=1,2,?,n-1) (1) 其中aj,bj,cj,dj 待定,并要使它满足: S(xj)=yj, S(xj-0)=S(xj+0), (j=2,?,n-1) (2) S'(xj-0)=S'(xj+0), S"(xj-0)=S"(xj+0), (j=2,?,n-1) (3) 式(2)、(3)共给出n+3(n-2)=4n-6 个条件, 需要待定4(n-1) 个系数,因此要唯一确定三次插值函数,还要附加2 个边界条件。通常由实际问题对三次样条插值在端点的状态要求给 出。常用边界的条件有以下3 类。 第 1 类边界条件:给定端点处的一阶导数值,S'(x1)=y1',S'(xn) =yn'。 第 2 类边界条件:给定端点处的二阶导数值,S"(x1)=y1",S"(xn) =yn"。特殊情况y1"=yn"=0,称为自然边界条件。 第 3 类边界条件是周期性条件,如果y=f(x)是以b-a 为周期的函 数,于是S(x) 在端点处满足条件S'(x1+0)=S'(xn-0),S"(x+0) =S"(xn-0)。 下以第 1 边界条件为例,利用节点处二阶导数来表示三次样条插值 告: Matlab中插值函数汇总和使用说明收藏 命令1 interp1 功能一维数据插值(表格查找。该命令对数据点之间计算内插值。它找出一元函数f(x在中间点的数值。其中函数f(x由所给数据决定。x:原始数据点 Y:原始数据点 xi:插值点 Yi:插值点 格式 (1yi = interp1(x,Y,xi 返回插值向量yi,每一元素对应于参量xi,同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵,则按Y 的每列计算。yi 是阶数为length(xi*size(Y,2的输出矩阵。 (2yi = interp1(Y,xi 假定x=1:N,其中N 为向量Y 的长度,或者为矩阵Y 的行数。 (3yi = interp1(x,Y,xi,method 用指定的算法计算插值: ’nearest’:最近邻点插值,直接完成计算; ’linear’:线性插值(缺省方式,直接完成计算; ’spline’:三次样条函数插值。对于该方法,命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函 数。命令spline 用它们执行三次样条函数插值; ’pchip’:分段三次Hermite 插值。对于该方法,命令interp1 调用函数p chip,用于对向量x 与y 执行分段三次内插值。该方法保留单调性与数据的外形; ’cubic’:与’pchip’操作相同; ’v5cubic’:在MATLAB 5.0 中的三次插值。 对于超出x 范围的xi 的分量,使用方法’nearest’、’linear’、’v5cubic’的插值算法,相应地将返回NaN。对其他的方法,interp1 将对超出的分量执行外插值算法。 (4yi = interp1(x,Y,xi,method,'extrap' 对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。 (5yi = interp1(x,Y,xi,method,extrapval 确定超出x 范围的xi 中的分量的外插值extrapval,其值通常取NaN 或0。 例1 1.>>x = 0:10; y = x.*sin(x; 2.>>xx = 0:.25:10; yy = interp1(x,y,xx; 3.>>plot(x,y,'kd',xx,yy 复制代码 例2 1.>> year = 1900:10:2010; 实验报告:牛顿差值多项式&三次样条 问题:在区间[-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)。 例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)); 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 学习报告—— 三次样条函数插值问题的讨论 班级:数学二班 学号: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、第二类边界条件:给定函数在端点处的二阶导数,即 数值计算方法作业 实验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 例:已知一组数据点,编写一程序求解三次样条插值函数满足 并针对下面一组具体实验数据 0.25 0.3 0.39 0.45 0.53 0.5000 0.5477 0.6245 0.6708 0.7280 求解,其中边界条件为. 1)三次样条插值自然边界条件源程序: function s=spline3(x,y,dy1,dyn) %x为节点,y为节点函数值,dy1,dyn分别为x=0.25,0.53处的二阶导 m=length(x);n=length(y); if m~=n error('x or y输入有误') return end h=zeros(1,n-1); h(n-1)=x(n)-x(n-1); for k=1:n-2 h(k)=x(k+1)-x(k); v(k)=h(k+1)/(h(k+1)+h(k)); u(k)=1-v(k); end g(1)=3*(y(2)-y(1))/h(1)-h(1)/2*dy1; g(n)=3*(y(n)-y(n-1))/h(n-1)+h(n-1)/2*dyn; for i=2:n-1 g(i)=3*(u(i-1)*(y(i+1)-y(i))/h(i)+v(i-1)*(y(i)-y(i-1))/h(i-1)); end for i=2:n-1; A(i,i-1)=v(i-1); A(i,i+1)=u(i-1); end A(n,n-1)=1; A(1,2)=1; A=A+2*eye(n); M=zhuigf(A,g); %调用函数,追赶法求M fprintf('三次样条(三对角)插值的函数表达式\n'); syms X; for k=1:n-1 fprintf('S%d--%d:\n',k,k+1); s(k)=(h(k)+2*(X-x(k)))./h(k).^3.*(X-x(k+1)).^2.*y(k)... +(h(k)-2*(X-x(k+1)))./h(k).^3.*(X-x(k)).^2.*y(k+1)... +(X-x(k)).*(X-x(k+1)).^2./h(k).^2*M(k)+(X-x(k+1)).*... (X-x(k)).^2./h(k).^2*M(k+1); end s=s.'; s=vpa(s,4); %画三次样条插值函数图像 for i=1:n-1 X=x(i):0.01:x(i+1); st=(h(i)+2*(X-x(i)))./(h(i)^3).*(X-x(i+1)).^2.*y(i)... +(h(i)-2.*(X-x(i+1)))./(h(i)^3).*(X-x(i)).^2.*y(i+1)... +(X-x(i)).*(X-x(i+1)).^2./h(i)^2*M(i)+(X-x(i+1)).*... (X-x(i)).^2./h(i)^2*M(i+1); plot(x,y,'o',X,st); hold on End plot(x,y); grid on %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %调用的函数: %追赶法 function M=zhuigf(A,g) n=length(A); L=eye(n); U=zeros(n); for i=1:n-1 U(i,i+1)=A(i,i+1); end U(1,1)=A(1,1); for i=2:n L(i,i-1)=A(i,i-1)/U(i-1,i-1); U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i); end Y(1)=g(1); for i=2:n Y(i)=g(i)-L(i,i-1)*Y(i-1); end M(n)=Y(n)/U(n,n); for i=n-1:-1:1 M(i)=(Y(i)-A(i,i+1)*M(i+1))/U(i,i); 告: Matlab中插值函数汇总和使用说明收藏 命令1 interp1 功能一维数据插值(表格查找)。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。x:原始数据点 Y:原始数据点 xi:插值点 Yi:插值点 格式 (1)yi = interp1(x,Y,xi) 返回插值向量yi,每一元素对应于参量xi,同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵,则按Y 的每列计算。yi 是阶数为length(xi)*size(Y,2)的输出矩阵。 (2)yi = interp1(Y,xi) 假定x=1:N,其中N 为向量Y 的长度,或者为矩阵Y 的行数。 (3)yi = interp1(x,Y,xi,method) 用指定的算法计算插值: ’nearest’:最近邻点插值,直接完成计算; ’linear’:线性插值(缺省方式),直接完成计算; ’spline’:三次样条函数插值。对于该方法,命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函 数。命令spline 用它们执行三次样条函数插值; ’pchip’:分段三次Hermite 插值。对于该方法,命令interp1 调用函数p chip,用于对向量x 与y 执行分段三次内插值。该方法保留单调性与数据的外形; ’cubic’:与’pchip’操作相同; ’v5cubic’:在MATLAB 5.0 中的三次插值。 对于超出x 范围的xi 的分量,使用方法’nearest’、’linear’、’v5cubic’的插值算法,相应地将返回NaN。对其他的方法,interp1 将对超出的分量执行外插值算法。 (4)yi = interp1(x,Y,xi,method,'extrap') 对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。 (5)yi = interp1(x,Y,xi,method,extrapval) 确定超出x 范围的xi 中的分量的外插值extrapval,其值通常取NaN 或0。 例1 1.>>x = 0:10; y = x.*sin(x); 2.>>xx = 0:.25:10; yy = interp1(x,y,xx); 3.>>plot(x,y,'kd',xx,yy) 复制代码 例2 1.>> year = 1900:10:2010; 2.>> product = [75.995 91.972 105.711 12 3.203 131.669 150.697 179.323 203.212 226.505 x0=[0 0.9211 1.8431 2.9497 3.8714 4.9781 5.9 7.0064 7.9286 8.9678 10.9542 12.0328 12.9544 13.8758 14.9822 15.9039 16.8261 17.9317 19.0375 19.9594 20.8392 22.9581 23.88 24.9869 25.9083]; >> >> y0=[14405 11180 10063 11012 8797 9992 8124 10160 8488 11018 19469 20196 18941 15903 18055 15646 13741 14962 16653 14496 14648 15225 15264 13708 9633]; >> x=0:0.1:25.9; >> y1=interp1(x0,y0,x,'spline'); >> pp1=csape(x0,y0); %样条插值工具箱函数 y2=ppval(pp1,x); %计算x对应的y值 pp2=csape(x0,y0,'second'); y3=ppval(pp2,x); xydata=[x',y1',y2',y3'] subplot(1,2,1) plot(x0,y0,'+',x,y1) title('Spline1') subplot(1,2,2) plot(x0,y0,'+',x,y2) title('Spline2') dx=diff(x); dy=diff(y2); dy_dx=dy./dx; dy_dx0=dy_dx(1) ytemp=y2(13<=x&x<=15); ymin=min(ytemp); xmin=x(y2==ymin); xymin_1315=[xmin,ymin] 实验报告:牛顿差值多项式&三次样条 问题:在区间[-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)。 《数值分析》课程设计 三次样条插值算法 院(系)名称信息工程学院 专业班级09普本信计1班 学号 学生姓名 指导教师 数值分析课程设计评阅书 课程设计任务书 2008—2009学年第二学期 专业班级:09普本信计1班学号:姓名: 课程设计名称:数值分析 设计题目:三次样条插值 完成期限:自2012 年 6 月8 日至2012 年 6 月13 日共 1 周 设计依据、要求及主要内容: 一、设计目的 熟练掌握三次样条插值算法的原理和推导过程,并且能够应用Matlab软件编写相应的程序和使用Matlab软件函数库软件。 二、设计要求 (1)用Matlab函数库中相应函数对选定的问题,求出具有一定精度的结果。 (2)使用所用的方法编写Matlab程序求解,对数值结果进行分析。 (3)对于使用多个方法解同一问题的,在界面上设计成菜单形式。 三、设计内容 首先构造三次样条插值函数的定义和一般特征,并对实例问题进行实例分析,并总结 四、参考文献 [1] 黄明游,冯果忱.数值分析[M].北京:高等教育出版社,2008. [2]马东升,雷勇军.数值计算方法[M].北京:机械工业出版社,2006. [3] 石博强,赵金.MATLAB数学计算与工程分析范例教程[M].北京:中国铁道出版社.2005. [4]郝红伟,MATLAB 6,北京,中国电力出版社,2001 [5]姜健飞,胡良剑,数值分析及其MATLAB实验,科学出版社,2004 [6]薛毅,数值分析实验,北京工业大学出版社,2005 计划答辩时间:2012年6月18日 指导教师(签字):教研室主任(签字):批准日期:年月 《数值分析课程设计》 报告 专业: 学号: 学生姓名: 指导教师: 7.掌握三次样条插值函数的构造方法,体会三次样条插值函数对被逼近函数的近似。 三次样条插值函数边界条件由实际问题对三次样条插值在端点的状态要求给出。以第1 边界条件为例, 用节点处二阶导数表示三次样条插值函数,用追赶法求解相关方程组。通过Matlab 编制三次样条函数的通用程序,可直接显示各区间段三次样条函数体表达式,计算出已给点插值并显示各区间分段曲线图。 引言 分段低次样条插值虽然计算简单、稳定性好、收敛性有保证且易在电子计算机上实现,但只能保证各小段曲线在连接处的连续性,不能保证整件曲线的光滑性。利用样条插值,既可保持分段低次插值多项式,又可提高插值函数光滑性。故给出分段三次样条插值的构造过程算法步骤,利用Matlab软件编写三次样条插值函数通用程序,并通过数值算例证明程序的正确性。 三次样条函数的定义及特征 定义:设[a,b] 上有插值节点,a=x1<x2<…x n=b,对应函数值为y1,y2,?y n。若函数S(x) 满足S(x j) =y j (j =1,2, ?,n ),S(x) 在[x j,x j+1] (j =1,2,?,n-1)上都是不高于三的多项式(为了与其对应j 从1 开始,在Matlab 中元素脚标从1 开始)。当S(x) 在[a,b] 具有二阶连续导数。则称S(x) 为三次样条插值函数。要求S(x) 只需在每个子区间[x j,x j+1] 上确定1 个三次多项式,设为: Sj(x)=ajx3+bjx2+cjx+dj, (j=1,2,?,n-1) (1) 其中a j,b j,c j,d j 待定,并要使它满足: S(x j)=y j, S(x j-0)=S(x j+0), (j=2,?,n-1) (2) S'(x j-0)=S'(x j+0), S"(x j-0)=S"(x j+0), (j=2,?,n-1) (3) 式(2)、(3)共给出n+3(n-2)=4n-6 个条件, 需要待定4(n-1) 个系数,因此要唯一确定三次插值函数,还要附加2 个边界条件。通常由实际问题对三次样条插值在端点的状态要求给出。常用边界的条件有以下3 类。 第1 类边界条件:给定端点处的一阶导数值,S'(x1)=y1',S'(x n)=y n'。 第 2 类边界条件:给定端点处的二阶导数值,S"(x1)=y1",S"(x n)=y n"。特殊情况y1"=y n"=0,称为自然边界条件。 第3 类边界条件是周期性条件,如果y=f(x)是以b-a 为周期的函数,于是S(x) 在端点处满足条件S'(x1+0)=S'(x n-0),S"(x+0)=S"(x n-0)。 下以第1 边界条件为例,利用节点处二阶导数来表示三次样条插值函数,给出具体的推导过程。 2 三次样条插值函数的推导过程 注意到S(x) 在[x j, x j+1](j=1,2,?,n-1)上是三次多项式,于是S"(x)在[x j, x j+1] 上是一次多项式,如果S"(x) 在[x j,x j+1](j=1,2,?,n-1)两端点上的值已知,设S"(x j)=M j,S"(x j+1)=M j+1,其中h j =x j+1-x j,对S"(x) 进行两次积分,则得到1 个具有2个任意常数A j,B j 的S(x) 表达式。对S"(x) 求两次积分 一.介绍 二.程序框图 开始 输入未知数X及 (xi,yi),i=0,1,…,n 计算步长H[i] 计算λ、μ、d 根据边界条件,求 解相应的方程得到 M1,…, Mn 将M代入原方程, 得到分段函数 结束 三.源码 syms h n=9;%插入节点数,可以根据题目更改 h=2/(n+1); u=0.5; v=0.5; f=inline('1/(1+25*x.^2)');%输入函数,这个也可以根据题目更改g=inline('3/h*((c-b)/h-(b-a)/h)','a','b','c','h'); for i=1:n+2 x(1)=-1; x(i+1)=x(i)+2/(n+1); y(i)=f(x(i)); end for i=1:n d(i)=g(y(i),y(i+1),y(i+2),h); end A=zeros(n,n); for i=1:n A(i,i)=2; end for i=1:n-1 A(i,i+1)=u; A(i+1,i)=v; end B=zeros(n,1); for i=1:n B(i,1)=d(i) end C=inv(A)*B for i=1:n M(i)=C(i,1); end x=(-1:h/50:1); k=1./(1+25*x.^2); cs=spline(x,k); plot(x,k,'r.'); hold on; ezplot('1/(1+25*x^2)',[-1 1]); title('三次样条插值曲线和f(x)曲线') 四.结果 系数矩阵 弯矩M 分段函数不同次幂x对应的系数 三次样条插值函数与原函数 三次样条插值多项式 ——计算物理实验作业四 陈万物理学2013级 主程序: clear,clc; format rat x = [1,4,9,16,25,36,49,64]; y = [1,2,3,4,5,6,7,8]; f1 = ; fn = 1/16; [a,b,c,d,M,S] = spline(x,y,f1,fn); 子程序1: function [a,b,c,d,M,S]=spline(x,y,f1,fn) % 三次样条插值函数 % x是插值节点的横坐标 % y是插值节点的纵坐标 % u是插值点的横坐标 % f1是左端点的一阶导数 % fn是右端点的一阶导数 % a是三对角矩阵对角线下边一行 % b是三对角矩阵对角线 % c是三对角矩阵对角线上边一行 % S是插值点的纵坐标 n = length(x); h = zeros(1,n-1); deltay = zeros(1,n); miu = zeros(1,n-1); lamda = zeros(1,n-1); d = zeros(1,n-1); for j = 1:n-1 h(j) = x(j+1)-x(j); deltay(j) = y(j+1)-y(j); end % 得到h矩阵 for j = 2:n-1 sumh = h(j-1) + h(j); miu(j) = h(j-1) / sumh; lamda(j) = h(j) / sumh; d(j) = 6*( deltay(j)/h(j)-(deltay(j-1)/h(j-1)))/sumh; end % 根据第一类边界条件,作如下规定 lamda(1) = 1; d(1) = 6*(deltay(1)/h(1)-f1)/h(1); miu(1) = 1; d(n) = 6*(fn-deltay(n-1)/h(n-1))/h(n-1);样条插值函数与应用
对样条函数及其插值问题的一点认识
(精品)数值分析课程设计-三次样条插值
Matlab中插值函数汇总和使用说明.
数值分析实验报告-插值、三次样条Word版
三次样条插值课后题集
三次样条插值方法的应用
关于三次样条插值函数的学习报告(研究生)资料
数值分析作业-三次样条插值
三次样条插值自然边界条件
Matlab中插值函数汇总和使用说明
三次样条插值函数matlab程序绝不坑爹
数值分析实验报告-插值、三次样条(教育教学)
数值分析课程设计--三次样条插值
数值分析课程设计(三次样条插值)
三次样条插值函数
三次样条插值多项式matlab