(精品)数值分析课程设计-三次样条插值
数值分析实验报告-插值、三次样条
实验报告:牛顿差值多项式&三次样条问题:在区间[-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 allclcx1=-1:0.2:1;y1=1./(1+25.*x1.^2);n=length(x1);f=y1(:);for j=2:nfor i=n:-1:jf(i)=(f(i)-f(i-1))/(x1(i)-x1(i-j+1));endendsyms F x p ;F(1)=1;p(1)=y1(1);for i=2:nF(i)=F(i-1)*(x-x1(i-1));p(i)=f(i)*F(i);endsyms PP=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 onxlabel('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)。
计算方法三次样条插值课程设计
摘要本文细致的讲解了三次样条插值函数的产生及在实际中解决的问题,通过MATLAB的程序编写,可以将复杂的计算省去,直接的给出了三次样条插值的结果与实际结果间的误差,验证实际结果和理论值的一致性。
避免了求解方程中的不必要计算,使求解效率得到显著的提高。
关键词插值函数三次样条插值 MATLAB1 三次样条插值函数概论当插值节点很多时,插值多项式的次数就会很高,这不仅增大了计算量,还会影响结果的精确度.虽然可以采用上述分段插值,但是主要缺点就是个分段接头处不光滑,插值函数的导数不连续,因此想构造这样的插值:既能分段的低次插值,又能保证接头处的光滑,就产生了三次样条插值函数.1.1定义设函数()f x 市区间[a,b]上的二次连续可微函数,在区间[a,b]上给处一个划分。
设函数()f x 是区间[a,b]上的一个划分011...n n a x x x x b -∆=<<<<=如果函数()S x 满足条件(1)在每个小区间1[,]k k x x +(k=1,2,….,n )上()S x 是一个部超过m 次的多项式。
(2)在节点k x (k=1,2,….,n )处具有m-1阶的连续导数。
(3)()()(0,1,2,...)j j s x f x j n ==1.2三次样条差值函数的构造由于三次样条插值我、函数s(x)的插值节点处的二阶导数存在,因此令各节点处的二阶导数为'()(0,1,...,)k s x m k n == (1.01)根据样条插值函数的定义,三次样条插值函数是s(x)在每一个小区间)1....,1,0](,[]1-=+n k x x k k 上市不超过三次的多项式。
在每一个小区间)1....,1,0](,[]1-=+n k x x k k 上,其二阶导数为线性函数,即''1111()k kk k k k k kx x x x s x m m x x x x ++++--=+-- (1.02)对式(1.02)积分两次,则得到k k k kk k k k k b x x a h x x m h x x m x s +-+-++=++)(6)(6)()(3131 (1.03)其中k k k k k b a x x h ,,1-=+为任意常数。
三次样条插值方法的应用
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在矩阵运算上的优势。
数学数值分析三次样条插值PPT课件
2.8.1 三次样条函数
定义 给定区间[a,b]的一个划分 a=x0<x1<…<xn=b, yi=f (xi) (i=0,1,…,n),如果函数S(x)满足: (1) S(xi )=yi (i=0,1,…,n); (2) 在每个小区间[xi, xi+1] (i=0,1,...,n-1)上是次数不超
S上且( xS与)(x相)的(邻x表节达x点j式的1 )为2两[hh个jj3转2角( x有关x j,)]故y j称为三h转j=x角j+方1-x程j 。
(
x
x
j
)2[hj 2( hj3
x
x j1 )] y j1
(x
x j1 )2 ( x h2j
xj)
mj
(x
x j )2( x h2j
x j1 )
m j1
则方程组化为:
2 1 2 2 2
m1 g1 1 f0
m2
g2
n2 2 n2 mn2 gn2
n1 2 mn1 gn1 n1 fn
第10页/共40页
2、已知 S( x0 ) f0, S( xn ) fn
2m0
m1
3
f
[x0 ,
x1 ]
h0 2
f0
第18页/共40页
S(
x)
M
j
(
x j1 6hj
x)3
M
j1
(x
x 6hj
j
)3
(
y
j
M jh2j 6
)
x
j1 hj
x
(
y
j1
M
j1h2j 6
)
x
x hj
数值分析——样条函数及三次样条插值
S k ( x )是[ xk , xk + 1 ]上的(两点)三次样条插值多项式, 满足
Sk ( x j ) = y j
x → xk
k = 0,1,2, ⋯ , n − 1; j = k , k + 1
lim S k ( x ) = lim S k − 1 ( x ) + −
一、三次样条插值函数
定义1. 定义
a ≤ x0 , x1 ,⋯ , xn ≤ b为区间[ a , b ]的一个分割
如果函数S ( x )在区间[ a , b ]上满足条件 :
( 1) S ( x ), S ′( x ), S ′′( x )都在区间[ a , b ]上连续 ,即 S ( x ) ∈ C 2 [ a , b ]
f ( x j ) = y j , j = 0 ,1,⋯ , n 如果S ( x )是f ( x )的三次样条插值函数, 则其必满足
S ( x j ) = y j , j = 0 ,1,⋯ , n lim S ( x ) = S ( x j ) = y j , j = 1,⋯ , n − 1 x→ x xlim S ′( x ) = S ′( x j ) = m j , j = 1,⋯ , n − 1 →x lim S ′′( x ) = S ′′( x j ), j = 1 ,⋯ , n − 1
+ x → xk
Sk ( x j ) = y j
k
k = 0,1, ⋯ , n − 1; j = k , k + 1
− x → xk
k −1
k = 1 , 2 ,⋯ , n − 1 k = 1,2 ,⋯ , n − 1 ------(8) k = 1, 2 ,⋯ , n − 1
数值计算方法三次样条插值
4.4 三次样条插值
令
A1
j1 (u )
(1
2
u
x hj
j 1
)(
u
xj hj
)2
A2
j (u )
(1
2
u
x hj
j
)(
u
x hj
j 1
)2
B1
j1 (u )
(u
u x j 1 )(
xj hj
)2
B2
j (u )
(u
x
j )(
u
x hj
j
)2
分段三次Hermite插值算法
I2(x)
I
n
(
x
)
x ( x0 , x1)
x ( x1, x2 ) ...... x ( xn1, xn )
其I中 j xxj1 xxj j yj1xxj xxjj 11yj yj1(xxj1)(yj yj1)/(xj xj1)
缺点:I(x)连续,但不光滑,精度较低,仅在 hm 1jan{xhj xj xj1}足够小才能较好。的逼近
ss((xxn0
) )
f f
( x0 ) (xn )
m0 mn
计算方法分段线性_三次样条插值
计算方法分段线性_三次样条插值分段线性和三次样条插值是两种常用的插值方法,在数值分析和插值问题中广泛使用。
1.分段线性插值分段线性插值是一种简单直观的插值方法,将插值区间划分为若干个子区间,在每个子区间上用线性函数进行插值。
假设给定的插值节点有n+1 个,节点为 (x0, y0), (x1, y1), ..., (xn, yn),并且满足 x0 <x1 < ... < xn。
则对于任意 xx 使得 x 在 [xi, xi+1] 之间,可以通过线性插值得到其函数值 yy,即:yy = yi + (xx - xi) * (yi+1 - yi) / (xi+1 - xi)分段线性插值方法简单易懂,适用于一些较简单的插值问题。
但是由于插值函数在节点之间是线性的,可能不能准确地反映出数据的特征,因此不适用于一些需要高精度的插值问题。
三次样条插值是一种更复杂、更精确的插值方法,将插值区间划分为若干个子区间,在每个子区间上用三次多项式进行插值。
三次样条插值方法的基本思想是找到一组三次多项式,满足在每个子区间内插值点的函数值和一阶导数值相等,并且两个相邻多项式在节点处的二阶导数值也相等。
具体的求解步骤如下:(1) 假设有 n+1 个插值节点 (x0, y0), (x1, y1), ..., (xn, yn),构造 n 个三次多项式,即每个多项式在 [xi, xi+1] 之间插值。
(2) 对每个子区间内的多项式进行插值,设第 i 个子区间的多项式为 Si(x) = ai + bi(x-xi) + ci(x-xi)^2 + di(x-xi)^3、将插值节点的函数值和一阶导数值代入多项式中,可以得到 n 个线性方程,利用这 n 个线性方程可以求解出 n 个子区间的系数。
(3)由于n个子区间的多项式必须在节点处一阶导数值相等,因此再设立n-1个方程,利用这些方程可以求解出n-1个子区间的二阶导数值。
(4)将求解得到的系数和二阶导数值代入每个子区间的多项式中,得到完整的三次样条插值函数。
数值分析三次样条插值
若取等距节点 hi = h, i = 1,…, n –1
i
h h
h
1 2
i
1 i
1 2
di
6 2h
yi 1
2 yi h
yi 1
3 h3
( yi1
2 yi
yi1 )
i 1, 2,, n
例1. 对于给定的节点及函数值
k 0123 xk 1 2 4 5 f (xk ) 1 3 4 2 求满足自然边界条件S(x0 ) S(xn ) 0的三次样条 插值函数S(x),并求f (3)的近似值
Mi1
( x xi )2 2hi 1
yi1 hi 1
yi
hi 1 6
( M i 1
Mi )
于是
Si( xi )
hi 3
Mi
yi
yi1 hi
hi 6
M i 1
Si1( xi )
hi 1 3
Mi
yi1 hi 1
yi
hi 1 6
M i 1
解: 由M关系式
k
hk
hk hk 1
k
hk 1 hk hk 1
1 k
1
2 3
1
1 3
2
1 3
2
2 3
di
6
yi1 hi1
yi
yi
yi hi
1
hi hi1 6 f [ xi1, xi , xi1]
三次样条插值算法详解
三次样条插值算法要求数据点数量较多,且在某些情况下可能存在数值不稳定性,如数据 点过多或数据点分布不均等情况。此外,该算法对于离散数据点的拟合效果可能不如其他 插值方法。
对未来研究的展望
01
02
03
改进算法稳定性
针对数值不稳定性问题, 未来研究可以探索改进算 法的数值稳定性,提高算 法的鲁棒性。
3
数据转换
对数据进行必要的转换,如标准化、归一化等, 以适应算法需求。
构建插值函数
确定插值节点
根据数据点确定插值节点,确保插值函数在节点处连续且光滑。
构造插值多项式
根据节点和数据点,构造三次多项式作为插值函数。
确定边界条件
根据实际情况确定插值函数的边界条件,如周期性、对称性等。
求解插值函数
求解线性方程组
06
结论
三次样条插值算法总结
适用性
三次样条插值算法适用于各种连续、光滑、可微的分段函数插值问题,尤其在处理具有复 杂变化趋势的数据时表现出色。
优点
该算法能够保证插值函数在分段连接处连续且具有二阶导数,从而在插值过程中保持数据 的平滑性和连续性。此外,三次样条插值算法具有简单、易实现的特点,且计算效率较高 。
根据数据点的数量和分布,合理分段,确保 拟合的精度和连续性。
求解线性方程组
使用高效的方法求解线性方程组,如高斯消 元法或迭代法。
结果输出
输出拟合得到的插值函数,以及相关的误差 分析和图表。
03
三次样条插值算法步骤
数据准备
1 2
数据收集
收集需要插值的原始数据点,确保数据准确可靠。
数据清洗
对数据进行预处理,如去除异常值、缺失值处理 等。
_三次样条插值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 车门曲线。
数值分析实验报告-插值、三次样条
实验报告:牛顿差值多项式&三次样条问题:在区间[-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 allclcx1=-1:0.2:1;y1=1./(1+25.*x1.^2);n=length(x1);f=y1(:);for j=2:nfor i=n:-1:jf(i)=(f(i)-f(i-1))/(x1(i)-x1(i-j+1));endendsyms F x p;F(1)=1;p(1)=y1(1);for i=2:nF(i)=F(i-1)*(x-x1(i-1));p(i)=f(i)*F(i);endsyms PP=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 onxlabel('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-14*x^3-16.855*x^2-6.6594e-16*x+1.0并且这里也能得到该牛顿插值多项式的在[-1,1]上的图形,并和原函数进行对比(见Fig.1)。
三次样条插值求导法
三次样条插值求导法三次样条插值法是一种常用的数值分析方法,用于近似插值实现平滑曲线的拟合。
它的优点在于可以保持原始数据的特性,同时能够降低数据间的噪声干扰,使得插值的结果更加准确。
本文将介绍三次样条插值法的原理、算法以及应用方面的指导意义。
首先,我们需要了解三次样条插值法的基本原理。
三次样条插值法通过在相邻数据点之间构造三次多项式来近似拟合原始数据。
这些三次多项式满足一定的光滑性条件,使得插值结果的曲线平滑而连续。
在三次样条插值中,每个数据点都对应一个三次多项式,并且相邻多项式之间的导数和二阶导数必须相等,以保证曲线的平滑性。
接下来,我们将介绍三次样条插值法的算法步骤。
首先,我们需要确定每个数据点对应的三次多项式。
为了满足光滑性条件,我们需要计算每个数据点处的导数值。
这可以通过求解一个线性方程组来实现,其中方程的个数等于数据点的个数。
解得导数值之后,我们就可以得到每个数据点对应的三次多项式的系数。
然后,我们需要利用这些系数来计算在数据点之间的插值结果。
为了实现这一点,我们可以利用三次多项式的性质,通过给定的数据点和对应的三次多项式系数,来计算在两个相邻数据点之间的插值结果。
最后,我们需要通过合理的选择数据点以及插值节点的间距,来获得更加准确的三次样条插值结果。
一般来说,数据点的选择应尽量满足曲线的变化趋势,以反映原始数据的特点。
此外,插值节点的间距也需要经过合理的选择,以保证插值结果的准确性。
三次样条插值法在实际应用中有着广泛的意义和指导价值。
首先,它可以用于光滑曲线的拟合,将离散的数据点进行连续化处理,使得数据的绘图和分析更加方便。
其次,它可以用于数据的插值预测,通过已有的数据点来预测未知数据点的取值。
此外,三次样条插值法还可以在数字图像处理中用于图像的平滑和插值填充,从而改善图像的质量和美观度。
综上所述,三次样条插值法是一种有效的数值分析方法,可以用于实现平滑曲线的拟合和数据的插值预测。
通过了解其原理、算法以及应用方面的指导意义,我们可以更好地理解和应用这一方法,从而提高数据处理和分析的准确性和效率。
数值分析三次样条插值函数
数值分析三次样条插值函数【问题】对函数f x =ex, x∈[0,1]构造等距节点的三次样条插值函数,对以下两种类型的样条函数1. 三次自然样条2. 满足S′ 0 =1,S′ 1 =e的样条并计算如下误差:max{ f x1 −S x1 ,i=1,…,N} i−i−i这里xi−1为每个小区间的中点。
对N=10,20,40比较以上两组节点的结果。
讨论你的结果。
【三次样条插值】在每一个区间[t1,t2],…,[tn−1,tn]上,S都是不同的三次多项式,我们把在[ti−1,ti]上表示S的多项式记为Si,从而,S0 x x∈[t0,t1]∈[t1,t2] S x = S1 x x…Sn−1 x x∈[tn−1,tn]通过在节点处函数值、一阶导数和二阶导数的连续性可以得到:Si−1 ti = yi= Si ti 1≤i≤ n−1Si−1′ ti = Si′ tix→ti+limS′′ x =zi=limS′′(x) x→ti−再给定z0和zn 的值就构成了4n个条件,而三次样条插值函数共4n个系数,故可以通过这4n个条件求解三次样条函数的系数,从而求得该三次样条插值函数。
特别的,当z0=zn=0 时称为自然三次样条。
文本预览:一、自然三次样条插值【自然三次样条插值算法】1.由上面的分析可知,求解三次样条函数实际上就是求解一个矩阵:u 1h 1h1u2h2h2u3…v1 z1 v2 z2 z3=v3 … z…hn−2 n−2 vn−2 z vn−1 un−1 n−1ih3…hn−3un−2hn−26…其中hi=ti+1−ti,ui=2(hi+hi−1),ui=h(yi+1−yi),vi=bi−bi−1 所以自然三层次样条插值的算法就是在得到端点的函数值,一次导数值和二次导数值,然后根据上述求解矩阵得到v,代入自然三次样条的表达式即可。
2.根据题目中所给出的误差估计,计算在区间中点处的最大误差。
【实验】通过Mathematica编写程序得到如下结果:N=101. 计算得到zi的值为:由此可以得到各个区间的自然三次样条插值函数。
三次样条插值法
鉴于高次插值不收敛又不稳定的特点,低次插值既具有收敛 性又具有稳定性,因此低次值更具有实用价值,但是低次插 值的光滑性较差,比如分段线性插值多项式在插值区间中仅 具有连续性,在插值节点处有棱角,一阶导数不存在;分段 三次Hermite插值多项式在插值区间中仅具有一阶导数即一 阶光滑性但不具备二阶光滑性,不能满足某些实际应用比如 汽车、轮船、飞机等的外形中流线形设计。样条是在二十世 纪初期经常用于图样设计的一种富有弹性的细长条,多个样 条互相弯曲连接后沿其边缘画出的曲线就是三次样条曲线。 后来数学上对其进行了抽象,定义了m次样条函数,并成为 数值逼近的重要研究分枝,进一步扩大了样条函数的应用范 围。
立即可得下式:
21
n m1 n mn1 2mn gn
其中:
h0 hn1 n , n 1 n h0 hn1 h0 hn1
y1 y0 yn yn1 g n 3 n n h0 hn 1
联合基本方程得一个广义三对角或周期三对角方程组:
S ( xn )
6 4 2 ( yn yn1 ) mn 1 mn M n 2 hn 1 hn 1 hn1
19
稍加整理得
y1 y0 h0 2m0 m1 3 M 0 g0 h0 2 yn yn 1 hn1 mn1 2mn 3 M n gn hn 1 2
hi i hi hi 1 hi 1 i 1 i hi hi 1
称为三转 角法基本 方程组
yi yi 1 yi 1 yi gi 3(i i ) hi 1 hi
共个n 1个方程, n 1个未知量
以上推导还没有考虑边界 条件!针对不同类型的三 次样条问题,就可以导出 不同的方程组!
数值分析课程设计--三次样条插值
《数值分析》课程设计三次样条插值算法院(系)名称信息工程学院专业班级 09普本信计1班学号 090111073学生姓名宣章然指导教师孔繁民2012年06月08日数值分析课程设计评阅书课程设计任务书2008—2009学年第二学期专业班级: 09普本信计1班学号: 060111060 姓名:宣章然课程设计名称:数值分析设计题目:三次样条插值完成期限:自 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日指导教师(签字):教研室主任(签字):批准日期:年月三次样条插值摘 要分段低次样条插值虽然计算简单、稳定性好、收敛性有保证且易在电子计算机上实现,但只能保证各小段曲线在连接处的连续性,不能保证整件曲线的光滑性。
利用样条插值,既可保持分段低次插值多项式,又可提高插值函数光滑性。
三次样条插值实验报告
x1=0:.01:1;y1=polyval(S1(1,:),x1-X(1)); x2=1:.01:2;y2=polyval(S1(2,:),x2-X(2)); x3=2:.01:3;y3=polyval(S1(3,:),x3-X(3)); x4=3:.01:4;y4=polyval(S1(4,:),x4-X(4)); x5=4:.01:5;y5=polyval(S1(5,:),x5-X(5)); x6=5:.01:6;y6=polyval(S1(6,:),x6-X(6)); >> plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,X,Y,'.') >> hold on >> x1=0:.01:1;y1=polyval(S2(1,:),x1-X(1)); x2=1:.01:2;y2=polyval(S2(2,:),x2-X(2)); x3=2:.01:3;y3=polyval(S2(3,:),x3-X(3)); x4=3:.01:4;y4=polyval(S2(4,:),x4-X(4)); x5=4:.01:5;y5=polyval(S2(5,:),x5-X(5)); x6=5:.01:6;y6=polyval(S2(6,:),x6-X(6)); >> plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,X,Y,'.') >> hold on >> x1=0:.01:1;y1=polyval(S3(1,:),x1-X(1)); x2=1:.01:2;y2=polyval(S3(2,:),x2-X(2)); x3=2:.01:3;y3=polyval(S3(3,:),x3-X(3)); x4=3:.01:4;y4=polyval(S3(4,:),x4-X(4)); x5=4:.01:5;y5=polyval(S3(5,:),x5-X(5)); x6=5:.01:6;y6=polyval(S3(6,:),x6-X(6)); >> plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,X,Y,'.')
数值分析作业-三次样条插值
数值计算方法作业实验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 边界条件为例,用节点处二阶导数表示三次样条插值函数,用追赶法求解相关方程组。
通过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 边界条件为例,利用节点处二阶导数来表示三次样条插值
函数,给出具体的推导过程。
2 三次样条插值函数的推导过程
注意到 S(x) 在 [xj, xj+1](j=1,2,⋯,n-1)上是三次多项式,于是S"(x)在[xj, xj+1] 上是一次多项式,如果S"(x) 在 [xj,xj+1](j=1,2,⋯,n-1)两端点上的值已知,设S"(xj)=Mj,S"(xj+1)=Mj+1,其中hj =xj+1-xj,对S"(x) 进行两次积分,则得到1 个具有2个任意常数Aj,Bj 的S(x) 表达式。
对S"(x) 求两次积分
4 三次样条插值函数的Matlab 程序设计
在 Matlab[3]环境下根据上述算法步骤进行编程,源程序如下:
function []=spline3(X,Y,dY,x0,m)
N=size(X,2);
s0=dY(1); sN=dY(2);
interval=0.025;
disp('x0为插值点')
x0
h=zeros(1,N-1);
for i=1:N-1 h(1,i)=X(i+1)-X(i); end
d(1,1)=6*((Y(1,2)-Y(1,1))/h(1,1)-s0)/h(1,1);
d(N,1)=6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1);
for i=2:N-1
d(i,1)=6*((Y(1,i+1)-Y(1,i))/h(1,i)-(Y(1,i)-Y(1,i-1))
/h(1,i-1))/ (h(1,i)+h(1,i-1)); endmu=zeros(1,N-1);
md=zeros(1,N-1);
md(1,N-1)=1; mu(1,1)=1;
for i=1:N-2
u=h(1,i+1)/(h(1,i)+h(1,i+1)); mu(1,i+1)=u;
md(1,i)=1-u; end
p(1,1)=2; q(1,1)=mu(1,1)/2;
for i=2:N-1
p(1,i)=2-md(1,i-1)*q(1,i-1); q(1,i)=mu(1,i)/p(1,i); end p(1,N)=2-md(1,N-1)*q(1,N-1);
y=zeros(1,N); y(1,1)=d(1)/2;
for i=2:N y(1,i)=(d(i)-md(1,i-1)*y(1,i-1))/p(1,i); end
x=zeros(1,N); x(1,N)=y(1,N);
for i=N-1:-1:1 x(1,i)=y(1,i)-q(1,i)*x(1,i+1); end
fprintf ('M为三对角方程的解\n'); M=x;
fprintf ('\n');
syms t;
digits (m);
for i=1:N-1
pp(i)=M(i)*(X(i+1)-t)^3/(6*h(i))+M(i+1)*(t-X(i))^3
/(6*h(i))+(Y(i)-M(i)*h(i)^2/6)*(X(i+1)-t)/h(i)+
(Y(i+1)-M(i+1)*h(i)^2/6)*(t-X(i))/h(i);
pp(i)=simplify(pp(i)); coeff=sym2poly(pp(i));
if length(coeff)~=4
tt=coeff(1:3); coeff(1:4)=0; coeff(2:4)=tt; end
if x0>X(i)&x0<X(i+1) L=i;
y0=coeff(1)*x0^3+coeff(2)*x0^2+coeff(3)*x0+coeff(4);
end
val=X(i):interval:X(i+1);
for k=1:length(val)
fval(k)=coeff(1)*val(k)^3+coeff(2)*val(k)^2+coeff(3)*val(k) +coeff(4); end
if mod(i,2)==1 plot(val,fval,'r+')
else plot(val,fval,'b.') end
hold on
clear val fval
ans=sym(coeff,'d');
ans=poly2sym(ans,'t');
fprintf('在区间[%f,%f]内\n',X(i),X(i+1));
fprintf('三次样条函数S(%d)=',i);
pretty(ans); end
fprintf ('x0所在区间为[%f,%f]\n',X(L),X(L+1));
fprintf ('函数在插值点x0=%f的值为\n',x0);
y0
程序中:X,Y 为输入结点,dY 为两端点一阶导数矩阵,x0 为待求插值点,m 为有效数字位数。
应用实例[1]:已知函数y=f(x)的函数值
如表1。
x0 为插值点x0=1.5;M 为三对角方程的解M=-5.0000 4.0000 4.0000 16.0000 在区间[-1.5,0.0]
内,三次样条函数S(1) t3 2t2 −1;在区间[0.0,1.0]内,三次样条函数S(2) 2t2 −1;在区间[1.0,2.0]内,三次样条函数S(3) 2t3 −4t2 6t −3。
5 结论
Matlab 环境下编写求解三次样条插值的通用程序,可直接显示各区间段三次样条函数的具体表达式,计算出已给点的插值,最后显示各区间分段曲线图,为三次样条插值函数的应用提供简便方法。