文档库 最新最全的文档下载
当前位置:文档库 › 三次样条函数的自动求法(学院+专业+学号)

三次样条函数的自动求法(学院+专业+学号)

三次样条函数的自动求法(学院+专业+学号)
三次样条函数的自动求法(学院+专业+学号)

三次样条函数的自动求法

摘要:在MATLAB 的The Spline Toolbox 中,没有给出三次样条函数表达的求法,可在教学过程中,或在实际问题中,我们需要知道样条插值函数的分段表达式。在现行数值分析教材中,一般都是通过解方程确定三次样条插值函数的表达式,但这种方法的工作量很大。在本文中,我们用MATLAB 语言编制了三个程序,给出在三种边界条件下,三次样条插值函数表达式的自动求法。

关键词:三次样条函数;边界条件;插值

0 引言

分段低次插值多项式具有计算简单、收敛性有保证、数值稳定等优点,但它不能保证整条曲线的光滑性甚至连续性,从而不能满足一些工程技术的要求。从20 世纪60 年代开始,由航空、造船等工程设计的需要而发展起来的样条插值方法,既保留了分段低次插值多项式的各种优点,又保证了插值函数的光滑性,已在许多领域里得到越来越广泛的应用。在教学过程中,或在实际问题中,我们需要知道样条插值函数的分段表达式。可在MATLAB 的The Spline Toolbox 中,没有直接给出三次样条函数表达式的求法,在现行数值分析教材中,一般都是在给定条件下,通过解方程而确定三次样条插值函数的表达式,尽管在计算过程中可借助数学软件来完成,但这种方法的工作量仍然很大。本文中,利用数学软件MATLAB ,我们给出了三次样条插值函数表达式的自动求法,这样不但解决了上述问题,而且给出了用数学软件解决实际问题的一个范例。

1 计算方法

定义对于给定的函数值

),,1,0)((n k x f y k k ==

其中b x x x a n =<<<= 10,如果函数)(x S 满足条件:

(1))(x S 在每个子区间[k k x x ,1-](k=1,2,n , )上都是不高于三次的多项式;

(2))(x S 、)(x S '、)(x S ''在[a,b]上都连续;

(3)),2,1,0()(n k y x S k ==。

则称)(x S 为函数)(x f 关于节点n x x x ,,,10 的三次样条插值函数。

要求三次样条插值函数)(x S ,只需在每个子区间[k k x x ,1-]上确定一个三次多项式k k k k k d x c x b x a x S +++=23)(

)(x S k 共有4个系数,

确定它们需要4个条件,因此要完全确定)(x S 共需4n 个条件。由)(x S 所满足的条件(1)、(2)、(3),可确定4n-2个条件,还缺少两个条件。这两个条件通常由实际问题对三次样条插值函数在端点的状态要求给出,称之为边界条件,常用的边界条件有以下三类。

第一类边界条件:给定端点处的一阶导数值

)(y x S '=',n n y x S '=')( 第二类边界条件:给定端点处的二阶导数值

)(y x S ''='',n n y x S ''='')( 当00

=''y ,0=''n y 时,)(x S 在端点处不受力,呈自然状态,故称之为自然边界条件。

第三类边界条件是周期性条件:如果)(x f y =是以b-a 为周期的函数,)(x S 也应是以b-a 为周期的函数,于是)(x S 在端点处满足条件

)0()0(00-'=+'x S x S

)0()0(-''=+''n n x S x S

设),,1,0()(n k M x S k k =='',记),,2,1(1n k x x h k k k =-=-,则可用k M 表示)(x S k :

k

k k k k k k k k k k k k k k k k h x x h M y h x x h M y h x x M h x x M x S 122113131)6()6(6)(6)()(-------+--+-+-=),,2,1],,[(1n k x x x k k =∈- (1)

若记

????

?????---+=-=+=+=-++++++)(611111111k k k k k k K k k k K k k k k k k k h y y h y y h h g h h h h h h μλμ (2) )1,,2,1(-=n k

在第一类边界条件下,可得确定0M 、1M 、……、n M 的线性方程

??????????

????????????=????????????????????????????????????????????----n n n n n n g g g g M M M M 1101101111

......212............212λμλμ (3) 其中???

????--'='--=-)(6)(61010110n n n n n n h y y y h g y h y y h g 在第二类边界条件下,可得确定0M 、1M 、……、n M 的线性方程组

??????????

????????????''-''-=???????????????????????????????????????????

?--------n n n n n n n n n y g g g y g M M M M 1122011122112222

1......2

2............22λμμλμλμλ (4) 其中:00y M ''=,n n y M ''=。在第三类边界条件下,可得确定1M 、……、0

M M n =的线性方程组

??????????

????????????=????????????????????????????????????????????----n n n n n n n n g g g g M M M M 1211211122

11......22............22μλλμλμμλ (5)

其中

????

?????---+=-=+=+=-)(6111011111n n n n n n n n n n n h y y h y y h h g h h h h h h μλμ 2 主要结果

在MATLAB 的The Spline Toolbox 中,没有直接给出三次样条函数表达式的求法,只给出其向量表示形式。根据(1)~(5)式,我们用MATLAB 语言编制了三个小程序,给出了在常用的三种边界条件下的三次样条插值函数表达式的自动求法。

命题 三种边界条件下的三次样条插值函数可分别由以下程序自动求得。 源程序 1 function S=spf1(x,y,xx,m)

%1 此函数为三次样条插值函数,计算满足边界条件一时的插值函数(此时需要定义一个系统变量’xx ’,即 syms xx ),函数的系数输出形式为分数,如果系数要以小数形式输出,请给出精度m 。

%2 如果给出的变量xx 不是系统变量,而是一组数值,则计算满足边界条件一时在xx 处的插值。

%3 输入时请将边界条件(即插值函数在两端点的一阶导数)%放在向量y 的第一个分量和最后一个分量。

%4 如果不输入边界条件,默认插值函数在两端点的一阶导数均为零。 n=length(x);

If n<2,error("There should be at least two data points."),end If any(diff(x)<0),[x,ind]=sort(x);else,ind=1:n;end

x=x(:);dx=diff(x);

If all(dx)==0,error("The data abscissae should be distinct."),end

[yd yn]=size(y);%if Y happen to be a column % matrix,change it to the expected row matrix.

If yn==1,yn=yd;y=reshape(y,1,yn);yd=1;end

If yn==n

notaknot=1;

elseif yn==n+2

notaknot=0;endslopes=y(:,[1 n+2]).':y(:,[1 n+2])=[];

Else

Error('Abscissa and ordinate vector should be of the same length.') End

yi=y(:.ind) ;dd=ones(1,yd);

Dx=diff(x);divdif=diff(yi)/dx(:,dd);

a=zeros(1,n-1);

c=a;

c=zeros(1,n-1);

a(1:n-2)=dx(1:n-2)/(dx(1:n-2)+dx(2:n-1));

c(2:n-1)=dx(2:n-1)/(dx(1:n-2)+dx(2:n-1));

d(2:n-1)=6*(divdif(2:n-1)-divdif(1:n-2))/(dx(1:n-2)+dx(2:n-1));

a(n-1)=1;

c(1)=1;

If notaknot

d(1)=6*(yi(2)-yi(1))/(dx(1)^2);

d(n)=6*(yi(n-1)-yi(n))/(dx(n-1)^2);

Else

d(1)=6*((yi(2)-yi(1))/dx(1)-endslopes(1))/dx(1);

d(n)=6*(endslopes(2)-(yi(n)-yi(n-1))/dx(n-1))/dx(n-1);

End

b(1:n)=2;

d=tridi(a,b,c,d);

If isnumeric(xx)==1

m=length(xx);

pp=0;

For k=1:m

For i=1:n-1

If(xx(k)<=x(i+1))&(xx(k)>=x(i))

pp(k)=d(i)*(x(i+1)-xx(k))^3/(6*dx(i))+d(i+1)*(xx(k)-x(i))^3/(6*dx(i)) +(y(i)-d(i)*dx(i)^2/6)*(x(i+1)-xx(k))/dx(i)+(y(i+1)-d(i+1)*dx(i)^2/6) *(xx(k)-x(i))/dx(i);

End

End

S=pp

End

Elseif(isnumeric(xx)+1==1)&(nargin==3)

For i=1:n-1

Pp(i)=d(i)*(x(i+1)-xx)^3/(6*dx(i))+d(i+1)*(xx(k)-x(i))^3/(6*dx(i))+(y (i)-(d(i)^2)/6)*(x(i+1)-xx(k))/dx(i)+(y(i+1)-(d(i+1)*dx(i)^2)/6)*(xx( k)-x(i))/dx(i);

Pp(i)=simplify(pp(i));

Fprintf(‘in[%g,%g]\n’,x(i),x(i+1));

Fprintf(‘S(%d)=’,i);

Pretty(pp(i));

End

Else

Digits(m);

For i=1:n-1

Pp(i)=d(i)*(x(i+1)-xx)^3/(6*dx(i))+d(i+1)*(xx-x(i))^3/(6*dx(i))+(y(i) -(d(i)^2)/6)*(x(i+1)-xx)/dx(i)+(y(i+1)-(d(i+1)*dx(i)^2)/6)*(xx-x(i))/ dx(i);

Sym2poly(pp(i));

Ans=sym(ans,’d’);

Poly2sym(ans,xx);

Fprintf(‘in[%g,%g]\n’,x(i),x(i+1));

Fprintf(‘S(%d)=’,i);

Pretty(ans);

End

End

3 应用实例

例 给定函数值k x =1,2,4,5,k y =f(k x )=1,3,4,2,求满足边界条件'S (1)=0,'

S (5)=1的三次样条插值函数。

解:(1)在命令窗口输入

xx=[1 2 4 5];y=[0 1 3 4 2 1]; syms x;spf1(xx,y,x) 输出: In[1,2] S(1)=

19835-14914x+225635x -39370

x In[2,4] S(2)=11435+19170x x+22235x -3314

x In[4,5] S(3)=-12267-184314x x-21585x +317370x spf1(xx,y,x,5)

输出: In[1,2] S(1)=-1.32863x +7.31432

x -10.643x+5.6571

In[2,4] S(2)=-0.214293x +0.628572x +2.7286x-3.2571

In[4,5] S(3)=2.47143x -31.6002x +131.64x-175.14

4 结论

利用此程序求三次样条插值函数的表达式是非常简单的事情,且它给出了两种表达形式。此程序还可用来求具体的插值。这样不但解决了上述问题,而且给出了用数学软件解决实际问题的一个范例。这给教学和实际应用提供了简单、方便、实用的工具。

样条插值函数与应用

样条插值函数及应用

摘要 样条函数具有广泛的应用,是现代函数论的一个十分活跃的分支,是计算方法的主要基础和工具之一,由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要,人们便更多地应用这个工具,也更深刻的认识了它的本质。 在实际问题中所遇到许多函数往往很复杂,有些甚至是很难找到解析表达式的。根据函数已有的数据来计算函数在一些新的点处的函数值,就是插值法所需要解决的问题。 插值法是数值逼近的重要方法之一,它是根据给定的自变量值和函数值,求取未知函数的近似值。早在一千多年前,我国科学家就在研究历法时就用到了线性插值和二次插值。而在实际问题中,有许多插值函数的曲线要求具有较高的光滑性,在整个曲线中,曲线不但不能有拐点,而且曲率也不能有突变。因此,对于插值函数必须二次连续可微且不变号 ,这就需要用到三次样条插值。 关键词三次样条函数;插值法

目录 引言 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 ++++ +≥??=??

三次样条插值课后题集

例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));

三次样条求导函数

三次样条求导函数 function ppd=ppder(pp) % pp=spline(x0,y0) [breaks,coefs,l]=unmkpp(pp); % breaks 是节点横坐标,coefs是每一段上的三次多项式的系数,l 是分的段 数 for n=1:l coefs1(n,:)=polyder(coefs(n,:)); end ppd=mkpp(breaks,coefs1,1); %ppd 是一个类似pp 的结构体数据,计算这样的分段多项式的值可以用 ppval(ppd,x) 举例: >> pp=spline(linspace(0,2*pi,50),sin(linspace(0,2*pi,50))) pp = form: 'pp' breaks: [1x50 double] coefs: [49x4 double] pieces: 49 order: 4 dim: 1 [breaks,coefs,l]=unmkpp(pp) breaks = Columns 1 through 12 0 0.1282 0.2565 0.3847 0.5129 0.6411 0.7694 … ….. Columns 49 through 50 6.1550 6.2832 coefs = -0.1643 -0.0007 1.0000 0 -0.1643 -0.0639 0.9918 0.1279

-0.1581 -0.1270 0.9673 0.2537 ……. 0.1453 -0.2731 -0.8381 0.5455 0.1546 -0.2172 -0.9010 0.4339 …… -0.1503 0.2457 0.8713 -0.4907 -0.1581 0.1879 0.9269 -0.3753 -0.1643 0.1270 0.9673 -0.2537 -0.1643 0.0639 0.9918 -0.1279 l = 49 >> pp1=mkpp(breaks,coefs,1) pp1 = form: 'pp' breaks: [1x50 double] coefs: [49x4 double] pieces: 49 order: 4 dim: 1 >>ppd=ppder(pp) ppd = form: 'pp' breaks: [1x50 double] coefs: [49x3 double] pieces: 49 order: 3 dim: 1 >>ppval(pp,3) ans = 0.1411 >> sin(3) ans = 0.1411

三次样条函数的自动求法(学院+专业+学号)

三次样条函数的自动求法 摘要:在MATLAB 的The Spline Toolbox 中,没有给出三次样条函数表达的求法,可在教学过程中,或在实际问题中,我们需要知道样条插值函数的分段表达式。在现行数值分析教材中,一般都是通过解方程确定三次样条插值函数的表达式,但这种方法的工作量很大。在本文中,我们用MATLAB 语言编制了三个程序,给出在三种边界条件下,三次样条插值函数表达式的自动求法。 关键词:三次样条函数;边界条件;插值 0 引言 分段低次插值多项式具有计算简单、收敛性有保证、数值稳定等优点,但它不能保证整条曲线的光滑性甚至连续性,从而不能满足一些工程技术的要求。从20 世纪60 年代开始,由航空、造船等工程设计的需要而发展起来的样条插值方法,既保留了分段低次插值多项式的各种优点,又保证了插值函数的光滑性,已在许多领域里得到越来越广泛的应用。在教学过程中,或在实际问题中,我们需要知道样条插值函数的分段表达式。可在MATLAB 的The Spline Toolbox 中,没有直接给出三次样条函数表达式的求法,在现行数值分析教材中,一般都是在给定条件下,通过解方程而确定三次样条插值函数的表达式,尽管在计算过程中可借助数学软件来完成,但这种方法的工作量仍然很大。本文中,利用数学软件MATLAB ,我们给出了三次样条插值函数表达式的自动求法,这样不但解决了上述问题,而且给出了用数学软件解决实际问题的一个范例。 1 计算方法 定义对于给定的函数值 ),,1,0)((n k x f y k k == 其中b x x x a n =<<<= 10,如果函数)(x S 满足条件: (1))(x S 在每个子区间[k k x x ,1-](k=1,2,n , )上都是不高于三次的多项式; (2))(x S 、)(x S '、)(x S ''在[a,b]上都连续; (3)),2,1,0()(n k y x S k ==。 则称)(x S 为函数)(x f 关于节点n x x x ,,,10 的三次样条插值函数。 要求三次样条插值函数)(x S ,只需在每个子区间[k k x x ,1-]上确定一个三次多项式k k k k k d x c x b x a x S +++=23)( )(x S k 共有4个系数, 确定它们需要4个条件,因此要完全确定)(x S 共需4n 个条件。由)(x S 所满足的条件(1)、(2)、(3),可确定4n-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

试求三次样条插值S(X)

给定数据表如下: 试求三次样条插值S(X),并满足条件: i)S’(0.25)=1.0000, S’(0.53)-0.6868; ii) S”(0.25)= S”(0.53)=0; 解: 由给定数据知: h0 =0.3-0.25 - 0.05 , h 1=0.39-0.30-0.09 h 2=0.45-0.39-0.06, h 3=0.53-0.45-0.08 由μ i=h i/(h i1+h i), λ i= h i/(h i1+h i) 得: μ1= 5/14 ; λ 1= 9/14 μ2= 3/5 ; λ 2= 2/5 μ3= 3/7 ; λ 3=4/7 0.25 0.5000 ﹨ ﹨ 1.0000 ∕﹨ 0.25 0.5000 ∕ -0.9200-f[x 0,x 0, x 1 ] ﹨∕ 0.9540 ∕﹨ 0.30 0.5477 -0.7193-f[x 0,x 1,x 2 ] ﹨∕

0.8533 ∕﹨ 0.39 0.6245 -0.5440-f[x1,x2,x 3 ] ﹨∕ 0.7717 ∕﹨ 0.45 0.6708 -0.4050-f[x 2,x 3,x 4 ] ﹨∕ 0.7150 ∕﹨ 0.53 0.7280 -0.3525-f[x 3,x 4,x 5 ] ﹨∕ 0.6868 ∕ 0.53 0.7280 i)已知一节导数边界条件,弯矩方程组 ┌┐┌┐ │ 2 1 │┌M 0 ┐│-0.9200 ︳ ︳5/14 2 9/14 ︳︳M ︳︳-0.7193 ︳ 1 ︳3/5 2 2/5 ︳︳M 2 ︳_6 ︳-0.5440︳ ︳ 3/7 2 4/7 ︳︳M ︳︳-0.4050 ︳ 3

关于三次样条插值函数的学习报告(研究生)资料

学习报告—— 三次样条函数插值问题的讨论 班级:数学二班 学号: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、第二类边界条件:给定函数在端点处的二阶导数,即

三次样条函数

计算方法实验报告 1、实验题目 三次样条插函数。 2、实验内容 三次样条插值是建立在Hermite 插值的基础上的。Hermite 插值是在一个区间上的插值,而三次样条插则是建立多个区间上插值,构造一个具有二阶光滑度的曲线,在求出给定点上对应的函数。本实验就是建立一个能根据三次样条插值函数求根的程序。 3、算法思想 给定一个区间,并把它分成n 等份,并且给出了每个结点对就的横坐标和纵坐标。利用程序输出给定插值点对应的值。横坐标设为:X 0, X 1, X 2, X 3, …X n 纵坐标为Y 0, Y 1, Y 2, …Y n ,设插点为u 。则令h k =X k+1-X k ,λk =1-+k k k h h h , μk =11--+k k k h h h , g k =3(1 11--+-+-k k k k k k k k h y y h y y λμ), 其中k=1,2,…,n-1 再根据第一类边界条件则可以确定公式6.16,再根据6.17解出方程中的m 向量,最后代入公式6.8求解。 4、源程序清单 #include #define N 21/*最大结点个数减一*/ void sanCi() { /*定义过程数据变量*/ float x[N],y[N],h[N]; /*横纵坐标及区间长度*/ float rr[N],uu[N],gg[N]; /*计算m 用的中间数组rr 、uu 、gg 分别对应:λ、μ、g 数组*/

float aa[N],bb[N],tt[N]; /*矩阵分解时用到的中间变量aa、bb、tt分别对应:α、β数组以及A=LU时中间矩阵*/ float mm[N]; /*最后要用到的系数m*/ int n,k,kv,chose; /* n为实际结点个数,k为下标,kv为最后确定k的值*/ float s,u; /*最后计算u对应的值*/ printf("请输入区间段数:"); scanf("%d",&n); /*输入结点个数*/ /*输入所有横坐标:*/ printf("输入所有横坐标:"); for(k=0; k<=n; k++) scanf("%f",&x[k]); /*输入对应纵坐标:*/ printf("输入对应纵坐标:"); for(k=0; k<=n; k++) scanf("%f",&y[k]); for(k=0; k

三次样条插值函数matlab程序绝不坑爹

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.1 定义 样条曲线是一个分段定义的公式。给定n+1个数据点,共有n个区间,三次样条方程满足以下条件: a. 在每个分段区间(i = 0, 1, …, n-1,x递增),都是一个三次多项式。 b. 满足(i = 0, 1, …, n ) c. ,导数,二阶导数在[a, b]区间都是连续的,即曲线是光滑的。 所以n个三次多项式分段可以写作: ,i = 0, 1, …, n-1 其中ai, bi, ci, di代表4n个未知系数。 1.2 求解 已知: a. n+1个数据点[xi, yi], i = 0, 1, …, n b. 每一分段都是三次多项式函数曲线 c. 节点达到二阶连续 d. 左右两端点处特性(自然边界,固定边界,非节点边界) 根据定点,求出每段样条曲线方程中的系数,即可得到每段曲线的具体表达式。 插值和连续性: , 其中i = 0, 1, …, n-1 微分连续性:

, 其中i = 0, 1, …, n-2 样条曲线的微分式: 将步长带入样条曲线的条件: a. 由(i = 0, 1, …, n-1)推出 b. 由(i = 0, 1, …, n-1)推出 c. 由(i = 0, 1, …, n-2)推出 由此可得: d. 由(i = 0, 1, …, n-2)推出 设,则 a. 可写为:

,推出 b. 将ci, di带入可得: c. 将bi, ci, di带入(i = 0, 1, …, n-2)可得: 端点条件 由i的取值范围可知,共有n-1个公式,但却有n+1个未知量m 。要想求解该方程组,还需另外两个式子。所以需要对两端点x0和xn的微分加些限制。选择不是唯一的,3种比较常用的限制如下。 a. 自由边界(Natural) 首尾两端没有受到任何让它们弯曲的力,即。具体表示为和 则要求解的方程组可写为: b. 固定边界(Clamped) 首尾两端点的微分值是被指定的,这里分别定为A和B。则可以推出

三次样条插值方法的应用

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 %dirivitive numbers given.

三次函数

三次函数 百科名片 三次函数 基本概念与性质形如y=ax^3+bx^2+cx+d(a≠0,b,c,d为常数)的函数叫做三次函数(cubics function)。三次函数的图像是一条曲线----回归式抛物线(不同于普通抛物线),具有比较特殊性。 目录 1二.零点求法1.盛金公式 12.盛金判别法 13.盛金定理 14.传统解法 三.三次函数性态的五个要点 1四.三次函数对称中心1.三次函数有对称中心 12.推广 五.其他性质 展开 编辑本段二.零点求法 求函数的零点可用盛金公式、盛金判别法、或传统解法盛金公式与盛金判别法及盛金定理的运用从这里向您介绍三次方程应用广泛。用根号解一元三次方程,虽然有著名的卡尔丹公式,并有相应的判别法,但使用卡尔丹公式解题比较复杂,缺乏直观性。范盛金推导出一套直接用a、b、c、d表达的较简明形式的一元三次方程的一般式新求根公式,并建立了新判别法。 1.盛金公式 一元三次方程aX^3+bX^2+cX+d=0,(a,b,c,d∈R,且a≠0)。重根判别式:A=b^2-3ac;B=bc-9ad;C=c^2-3bd,总判别式:

Δ=B^2-4AC。当A=B=0时,盛金公式①:X1=X2=X3=-b/(3a)=-c/b=-3d/c。当Δ=B^2-4AC>0时,盛金公式②:X1=(-b-(Y1)^(1/3)-(Y2)^(1/3))/(3a);X2,3=(-2b+(Y1)^(1/3)+(Y2)^(1/3))/(6a)±i3^(1/2)((Y1)^(1/3)-(Y2)^(1/3))/(6a),其中Y1,2=Ab+3a(-B±(B^2-4AC)^(1/2))/2,i^2=-1。当Δ=B^2-4AC=0时,盛金公式③:X1=-b/a+K;X2=X3=-K/2,其中K=B/A,(A≠0)。当Δ=B^2-4AC<0时,盛金公式④:X1= (-b-2A^(1/2)cos(θ/3))/(3a);X2,3= (-b+A^(1/2)(cos(θ/3)±3^(1/2)sin(θ/3)))/(3a),其中θ=arccosT,T= (2Ab-3aB)/(2A^(3/2)),(A>0,-10时,方程有一个实根和一对共轭虚根;③:当Δ=B^2-4AC=0时,方程有三个实根,其中有一个两重根;④:当Δ=B^2-4AC<0时,方程有三个不相等的实根。 3.盛金定理 当b=0,c=0时,盛金公式①无意义;当A=0时,盛金公式③无意义;当A≤0时,盛金公式④无意义;当T<-1或T>1时,盛金公式④无意义。当b=0,c=0时,盛金公式①是否成立?盛金公式③与盛金公式④是否存在A≤0的值?盛金公式④是否存在T<-1或T>1的值?盛金定理给出如下回答:盛金定理1:当A=B=0时,若b=0,则必定有c=d=0(此时,方程有一个三重实根0,盛金公式①仍成立)。盛金定理2:当A=B=0时,若b≠0,则必定有c≠0(此时,适用盛金公式①解题)。盛金定理3:当A=B=0时,则必定有C=0(此时,适用盛金公式①解题)。盛金定理4:当A=0时,若B≠0,则必定有Δ>0(此时,适用盛金公式②解题)。盛金定理5:当A<0时,则必定有Δ>0(此时,适用盛金公式②解题)。盛金定理6:当Δ=0时,若B=0,则必定有A=0(此时,适用盛金公式①解题)。盛金定理7:当Δ=0时,若B≠0,盛金公式③一定不存在A≤0的值(此时,适用盛金公式③解题)。盛金定理8:当Δ<0时,盛金公式④一定不存在A≤0的值。(此时,适用盛金公式④解题)。盛金定理9:当Δ<0时,盛金公式④一定不存在T≤-1或T≥1的值,即T出现的值必定是-10时,不一定有A<0。盛金定理表明:盛金公式始终保持有意义。任意实系数的一元三次方程都可以运用盛金公式直观求解。当Δ=0(d≠0)时,使用卡尔丹公式解题仍存在开立方。与卡尔丹公式相比较,盛金公式的表达形式较简明,使用盛金公式解题较直观、效率较高;盛金判别法判别方程的解较直观。重根判别式A=b^2-3ac;B=bc-9ad;C=c^2-3bd是最简明的式子,由A、B、C构成的总判别式Δ=B^2-4AC也是最简明的式子(是非常美妙的式子),其形状与一元二次方程的根的判别式相同;盛金公式②中的式子(-B±(B^2-4AC)^(1/2))/2具有一元二次方程求根公式的形式,这些表达形式体现了数学的有序、对称、和谐与简洁美。

三次样条插值多项式matlab

三次样条插值多项式 ——计算物理实验作业四 陈万物理学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);

三次样条曲线插补模拟

主程序M文件(King.m): function varargout = King(varargin) % KING MATLAB code for King.fig % KING, by itself, creates a new KING or raises the existing % singleton*. % % H = KING returns the handle to a new KING or the handle to % the existing singleton*. % % KING('CALLBACK',hObject,eventData,handles,...) calls the local % function named CALLBACK in KING.M with the given input arguments. % % KING('Property','Value',...) creates a new KING or raises the % existing singleton*. Starting from the left, property value pairs are % applied to the GUI before King_OpeningFcn gets called. An % unrecognized property name or invalid value makes property application % stop. All inputs are passed to King_OpeningFcn via varargin. % % *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one % instance to run (singleton)". % % See also: GUIDE, GUIDATA, GUIHANDLES % Edit the above text to modify the response to help King % Last Modified by GUIDE v2.5 06-May-2013 21:39:33 % Begin initialization code - DO NOT EDIT gui_Singleton = 1; gui_State = struct('gui_Name', mfilename, ... 'gui_Singleton', gui_Singleton, ... 'gui_OpeningFcn', @King_OpeningFcn, ... 'gui_OutputFcn', @King_OutputFcn, ... 'gui_LayoutFcn', [] , ... 'gui_Callback', []);

数值分析作业-三次样条插值教学提纲

数值分析作业-三次样 条插值

数值计算方法作业 实验4.3 三次样条差值函数 实验目的: 掌握三次样条插值函数的三弯矩方法。 实验函数: dt e x f x t ? ∞ -- =2 221)(π 求f(0.13)和f(0.36)的近似值 实验内容: (1) 编程实现求三次样条插值函数的算法,分别考虑不同的边界条件; (2) 计算各插值节点的弯矩值; (3) 在同一坐标系中绘制函数f(x),插值多项式,三次样条插值多项式的曲 线比较插值结果。 实验4.5 三次样条差值函数的收敛性 实验目的: 多项式插值不一定是收敛的,即插值的节点多,效果不一定好。对三次样条插值函数如何呢?理论上证明三次样条插值函数的收敛性是比较困难的,通过本实验可以证明这一理论结果。

实验内容: 按照一定的规则分别选择等距或非等距的插值节点,并不断增加插值节点的个数。 实验要求: (1) 随着节点个数的增加,比较被逼近函数和三样条插值函数的误差变化情 况,分析所得结果并与拉格朗日插值多项式比较; (2) 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子, 考虑如下例子:某汽车制造商根据三次样条插值函数设计车门曲线,其中一段数据如下: 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 k y ' 0.8 0.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

三次样条拟合典型实例

1设计目的、要求 对龙格函数2 2511 )(x x f += 在区间[-1,1]上取10=n 的等距节点,分别作多项式插值、三次样条插值和三次曲线拟合,画出)(x f 及各逼近函数的图形,比较各结果。 2设计原理 (1) 多项式插值:利用拉格朗日多项式插值的方法,其主要原理是拉格朗日多项 式,即: 01,,...,n x x x 表示待插值函数的1n +个节点, 0()()n n j k k j j k L x y l x y ===∑,其中0,1,...,j n =; 011011()...()()...() ()()...()...()...() k k n k k k k k k k n x x x x x x x x l x x x x x x x x x -+-+----= ---- (2) 三次样条插值:三次样条插值有三种方法,在本例中,我们选择第一边界条 件下的样条插值,即两端一阶导数已知的插值方法: 00'()'S x f = '()'n n S x f = (3)三次曲线拟合:本题中采用最小二乘法的三次多项式拟合。最小二乘拟合是 利用已知的数据得出一条直线或者曲线,使之在坐标系上与已知数据之间的距离的平方和最小。在本题中,n= 10,故有11个点,以这11个点的x 和 y 值为已知数据,进行三次多项式拟合,设该多项式为23432xi i i i p a a x a x ax =+++,该拟合曲线只需2[]xi i p y -∑的值最小即可。 3采用软件、设备 计算机、matlab 软件

4设计内容 1、多项式插值: 在区间[]1,1-上取10=n 的等距节点,带入拉格朗日插值多项式中,求出各个节点的插值,并利用matlab 软件建立m 函数,画出其图形。 在matlab 中建立一个lagrange.m 文件,里面代码如下: %lagrange 函数 function y=lagrange(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end 建立一个polynomial.m 文件,用于多项式插值的实现,代码如下: %lagrange 插值 x=[-1:0.2:1]; y=1./(1+25*x.^2); x0=[-1:0.02:1]; y0=lagrange(x,y,x0); y1=1./(1+25*x0.^2); plot(x0,y0,'--r') %插值曲线 hold on %原曲线 plot(x0,y1,'-b') 运行duoxiangshi.m 文件,得到如下图形:

相关文档