计算方法三次样条插值课程设计

合集下载

(精品)数值分析课程设计-三次样条插值

(精品)数值分析课程设计-三次样条插值

《数值分析课程设计-三次样条插值》报告掌握三次样条插值函数的构造方法,体会三次样条插值函数对被逼近函数的近似。

三次样条插值函数边界条件由实际问题对三次样条插值在端点的状态要求给出。

以第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个边界条件。

数值分析课程设计报告书三次样条插值的三弯矩法

数值分析课程设计报告书三次样条插值的三弯矩法

数值分析课程设计报告书院系名称:学生姓名:专业名称:班级:时间:实验一 三次样条插值的三弯矩法一、实验目的已知数据i x ,()i i y f x =,0,,i n =及边界条件()n j x y j j 1,0),(2=,求)(x f 的三次样条插值函数)(x S .要求输出用追赶法解出的弯矩向量0[,,]n M M M =及()(),0,,,0,1,2k i S t i m k ==的值.画出)(x S y =的图形,图形中描出插值点(,)i i x y 及(,())i i t S t 分别用‘o ’和‘*’标记.二、实验原理1.用追赶法求解第二类边界条件的三弯矩方程:0010012111121111[,,]21[,,]26[,,]212[,,]n n n n n n n n n n f x x x M f x x x M M f x x x M f x x x μλμλ------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 其中1111,,j jj j j j j j j j j h h h x x h h h h μλ-+--===-++.2.得出样条函数表达式:332211111()()()()()6666j j j j j j j j j j j j j j j jx x x x M h x x M h x x S x M M y y h h h h +++++----=++-+-. 3.计算(k)(),0,,,0,1,2i S t i m k ==.三、实验结果所用数据:x=[-2.223,-1.987,-1.8465,-1.292,-1.2266,-1.1056,-0.8662,-0.6594,-0.2671,-0.0452,0.5385,1.2564,1.4398,1.5415,1.7646,1.9678,2.236];y=[0.83995,1.1696,1.3141,1.6992,1.7312,1.7847,1.8708,1.9262,1.9881,1.9997,1.9511,1.7169,1.618,1.5543,1.3871,1.191,0.81662];d2s1= -4.5000;d2sn= -4.8967; %第二种边界条件t=[-2.223,-1.9443,-1.6656,-1.3869,-1.1083,-0.82956,-0.55088,-0.27219,0.0065,0.28519,0.56387,0.84256,1.1212,1.3999,1.6786,2.236]; ;(指定计算点)计算结果:-2.5-2-1.5-1-0.500.51 1.52 2.50.811.21.41.61.82四、实验分析通过实验结果我们,知道三弯矩法求出满足初始条件的三次样条函数,与其他插值函数的构造相比,三次样条插值法的计算量要小得多。

三次样条插值的方法和思路

三次样条插值的方法和思路

三次样条插值的方法和思路摘要:1.三次样条插值的基本概念2.三次样条插值的数学原理3.三次样条插值的实现步骤4.三次样条插值的优缺点5.三次样条插值在实际应用中的案例正文:在日常的科学研究和工程应用中,我们经常会遇到需要对一组数据进行插值的问题。

插值方法有很多,其中三次样条插值是一种常见且有效的方法。

本文将从基本概念、数学原理、实现步骤、优缺点以及实际应用案例等方面,全面介绍三次样条插值的方法和思路。

一、三次样条插值的基本概念三次样条插值(Cubic Spline Interpolation)是一种基于分段多项式的插值方法。

它通过在各个节点上构建一条三次多项式曲线,使得这条曲线在节点之间满足插值条件,从而达到拟合数据的目的。

二、三次样条插值的数学原理三次样条插值的数学原理可以分为两个部分:一是分段三次多项式的构建,二是插值条件的满足。

1.分段三次多项式的构建假设有一组数据点序列为(x0,y0),(x1,y1),(x2,y2),(x3,y3),我们可以将这些数据点连接起来,构建一条分段三次多项式曲线。

分段三次多项式在每个子区间上都是一个三次多项式,它们之间通过节点值进行连接。

2.插值条件的满足为了使分段三次多项式在节点之间满足插值条件,我们需要在每个子区间上满足以下四个条件:(1)端点条件:三次多项式在区间的端点上分别等于节点值;(2)二阶导数条件:三次多项式在区间内的二阶导数等于节点间的斜率;(3)三阶导数条件:三次多项式在区间内的三阶导数等于节点间的曲率;(4)内部点条件:三次多项式在区间内部满足插值函数的连续性。

通过求解这四个条件,我们可以得到分段三次多项式的系数,从而实现插值。

三、三次样条插值的实现步骤1.确定插值节点:根据数据点的位置,选取合适的节点;2.构建分段三次多项式:根据节点值和插值条件,求解分段三次多项式的系数;3.计算插值结果:将待插值点的横坐标代入分段三次多项式,得到插值结果。

样条函数及三次样条插值PPT课件

样条函数及三次样条插值PPT课件

(x)
lim
x xk
Sk 1( x)
lim
x
x
k
Sk (x)
lim
x
x
k
Sk1( x)
k 1,2,,n 1
------(4)
lim
x
x
k
Sk( x)
lim
x
x
k
Sk1( x)
共4n 2个条件
5
Sk (x)是[xk , xk 1 ]上的三次样条插值多项式,应有4个待定的系数 即要确定S(x)必须确定4n个待定的系数 少两个条件 并且我们不能只对插值函数在中间节点的状态进行限制 也要对插值多项式在两端点的状态加以要求 也就是所谓的边界条件:
例. 使用不同的插值方法于函数
y
1
1 x2
x [5,5]
最后,介绍一个有用的结论
定理 . 设f (x) C 2[a,b], S(x)是以xk (k 0,1,, n)
为节点, 满足任意边界条件的三次样条插值函数,
设hi
xi 1
xi
,
h
max
0in1
hi
,
min
0in1
hi
,
则当 h
c 时
S(x)和S(x)在[a,b]上一致收敛到f (x)和f (x)
------(6)
13
由(11)式,可知
S0( x0
)
6( x0
x1 h03
2 x0
) ( y1
y0 )
6 x0
2 x0 h02
4 x1
m0
6 x0
4 x0 h02
2 x1
m1
6 h02
(

三次样条插值的求解

三次样条插值的求解

三次样条插值的求解摘要:分段低次插值虽然解决了高次插值的振荡现象和数值不稳定现象,使得插值多项式具有一致收敛性,保证了插值函数整体的连续性,但在函数插值节点处不能很好地保证光滑性要求,这在某些要求光滑性的工程应用中是不能接受的。

如飞机的机翼一般要求使用流线形设计,以减少空气阻力,还有船体放样等的型值线,往往要求有二阶光滑度(即有二阶连续导数)。

因此,在分段插值的基础上,引进了一种新的插值方法,在保证原方法的收敛性和稳定性的同时,又使得函数具有较高的光滑性的样条插值。

关键字:三转角方程 三弯矩阵方程0. 引言1,三次样条函数定义1:若函数2()[,]S x a b C ∈,且在每个小区间上1,j j x x +⎡⎤⎦⎣上是三次多项式,其中01n a x x x b ⋯=<<<= 是给定节点,则称()s x 是节点01,,,n x x x ⋯上的三次样条函数。

若节点j x 上 给定函数值()j j y f x =(0,1,)j n ⋯= ,且()j j s x y = (0,1,)j n ⋯= (1.1)成立,则称 ()s x 为三次样条差值函数。

从定义知,要求出()s x ,在每个应小区间1[,]j j x x + 上确定4个待定系数,共有 n 个小区间,故应确定4n 个参数,根据()s x 在[,]a b 上二阶导数连续,在节点()1,2,3,,1j x j n ⋯=-处应满足连续性条件(0)(0),j j s x s x -=+ ''(0)(0),j j s x s x -=+''''(0)(0)j j s x s x -=+ (1.2) 共有 3n-3个条件,再加上()s x 满足插值条件(1.1),共有4n-2个条件,因此还需要2个条件才能确定()s x 。

通常可在区间[,]a b 端点0,n a x b x ==上各加一个条件(称边界条件),边界条件可根据实际的问题要求给定。

三次样条插值计算算法

三次样条插值计算算法

/* 三次样条插值计算算法*/#include "math.h "#include "stdio.h "#include "stdlib.h "/*N:已知节点数N+1R:欲求插值点数R+1x,y为给定函数f(x)的节点值{x(i)} (x(i) <x(i+1)) ,以及相应的函数值{f(i)} 0 <=i <=NP0=f(x0)的二阶导数;Pn=f(xn)的二阶导数u:存插值点{u(i)} 0 <=i <=R求得的结果s(ui)放入s[R+1] 0 <=i <=R返回0表示成功,1表示失败*/int SPL(int N,int R,double x[],double y[],double P0,double Pn,double u[],double s[]){/*声明局部变量*/double *h; /*存放步长:{hi} 0 <=i <=N-1 */double *a; /*存放系数矩阵{ai} 1 <=i <=N ;分量0没有利用*/ double *c; /*先存放系数矩阵{ci} 后存放{Bi} 0 <=i <=N-1 */double *g; /*先存放方程组右端项{gi} 后存放求解中间结果{yi} 0 <=i <=N */double *af; /*存放系数矩阵{a(f)i} 1 <=i <=N ;*/double *ba; /*存放中间结果0 <=i <=N-1*/double *m; /*存放方程组的解{m(i)} 0 <=i <=N ;*/int i,k;double p1,p2,p3,p4;/*分配空间*/if(!(h=(double*)malloc(N*sizeof(double)))) exit(1);if(!(a=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(c=(double*)malloc(N*sizeof(double)))) exit(1);if(!(g=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(af=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(ba=(double*)malloc((N)*sizeof(double)))) exit(1);if(!(m=(double*)malloc((N+1)*sizeof(double)))) exit(1);/*第一步:计算方程组的系数*/for(k=0;k <N;k++)h[k]=x[k+1]-x[k];for(k=1;k <N;k++)a[k]=h[k]/(h[k]+h[k-1]);for(k=1;k <N;k++)c[k]=1-a[k];for(k=1;k <N;k++)g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]); c[0]=a[N]=1;g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;/*第二步:用追赶法解方程组求{m(i)} */ba[0]=c[0]/2;g[0]=g[0]/2;for(i=1;i <N;i++){af[i]=2-a[i]*ba[i-1];g[i]=(g[i]-a[i]*g[i-1])/af[i];ba[i]=c[i]/af[i];}af[N]=2-a[N]*ba[N-1];g[N]=(g[N]-a[N]*g[N-1])/af[N];m[N]=g[N]; /*P110 公式:6.32*/ for(i=N-1;i> =0;i--)m[i]=g[i]-ba[i]*m[i+1];/*第三步:求值*/for(i=0;i <=R;i++){/*判断u(i)属于哪一个子区间,即确定k */if(u[i] <x[0] || u[i]> x[N]){/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m);return 1;}k=0;while(u[i]> x[k+1])k++;//p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3); //p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3);p1=(h[k]+2*(u[i]-x[k]))*pow((u[i]-x[k+1]),2)*y[k]/pow(h[k],3);p2=(h[k]-2*(u[i]-x[k+1]))*pow((u[i]-x[k]),2)*y[k+1]/pow(h[k],3); p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2);p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2);s[i]=p1+p2+p3+p4;}/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m);return 0;}void main(){int N,R;double *x,*y,*u,*s;double P0,Pn;int i;/*验证算法:*/N=7;R=6;/*分配空间*/if(!(x=(double*)malloc((N+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(y=(double*)malloc((N+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(u=(double*)malloc((R+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(s=(double*)malloc((R+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}x[0]=0.5;x[1]=0.7;x[2]=0.9;x[3]=1.1;x[4]=1.3;x[5]=1.5;x[6]=1.7;x[7]=1.9;y[0]=0.4794;y[1]=0.6442;y[2]=0.7833;y[3]=0.8912;y[4]=0.9636;y[5]=0.9975;y[6]=0.9917;y[7]=0.9 463;u[0]=0.6;u[1]=0.8;u[2]=1.0;u[3]=1.2;u[4]=1.4;u[5]=1.6;u[6]=1.8;P0=-0.4794;Pn=-0.9463;if(!SPL( N, R, x, y, P0, Pn, u, s)){/*打印结果*/printf( "\nx= ");for(i=0;i <=N;i++)printf( "%8.1f ",x[i]);printf( "\ny= ");for(i=0;i <=N;i++)printf( "%8.4f ",y[i]);printf( "\n\nu= ");for(i=0;i <=R;i++)printf( "%9.2f ",u[i]);printf( "\ns= ");for(i=0;i <=R;i++)printf( "%9.5f ",s[i]);printf( "\nsin= ");for(i=0;i <=R;i++)printf( "%9.5f ",sin(u[i]));}/*释放空间*/free(x);free(y);free(u);free(s);}/* 测试数据来自课本55页例5 《数值分析》清华大学出版社第四版*/ //输入327.7 4.128 4.329 4.130 3.013.0 -4.0//输出输出三次样条插值函数:1: [27.7 , 28]13.07*(x - 28)^3 + 0.22*(x - 27.7)^3+ 14.84*(28 - x) + 14.31*(x - 27.7)2: [28 , 29]0.066*(29 - x)^3 + 0.1383*(x - 28)^3+ 4.234*(29 - x) + 3.962*(x - 28)3: [29 , 30]0.1383*(30 - x)^3 - 1.519*(x - 29)^3+ 3.962*(30 - x) + 4.519*(x - 29)//三次样条插值函数#include<iostream>#include<iomanip>using namespace std;const int MAX = 50;float x[MAX], y[MAX], h[MAX];float c[MAX], a[MAX], fxym[MAX];float f(int x1, int x2, int x3){float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);return (a - b)/(x[x3] - x[x1]);} //求差分void cal_m(int n){ //用追赶法求解出弯矩向量M……float B[MAX];B[0] = c[0] / 2;for(int i = 1; i < n; i++)B[i] = c[i] / (2 - a[i]*B[i-1]);fxym[0] = fxym[0] / 2;for(i = 1; i <= n; i++)fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);for(i = n-1; i >= 0; i--)fxym[i] = fxym[i] - B[i]*fxym[i+1];}void printout(int n);int main(){int n,i; char ch;do{cout<<"Please put in the number of the dots:";cin>>n;for(i = 0; i <= n; i++){cout<<"Please put in X"<<i<<':';cin>>x[i]; //cout<<endl;cout<<"Please put in Y"<<i<<':';cin>>y[i]; //cout<<endl;}for(i = 0; i < n; i++) //求步长h[i] = x[i+1] - x[i];cout<<"Please 输入边界条件\n 1: 已知两端的一阶导数\n 2:两端的二阶导数已知\n 默认:自然边界条件\n";int t;float f0, f1;cin>>t;switch(t){case 1:cout<<"Please put in Y0\' Y"<<n<<"\'\n";cin>>f0>>f1;c[0] = 1; a[n] = 1;fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];break;case 2:cout<<"Please put in Y0\" Y"<<n<<"\"\n";cin>>f0>>f1;c[0] = a[n] = 0;fxym[0] = 2*f0; fxym[n] = 2*f1;break;default:cout<<"不可用\n";//待定};//switchfor(i = 1; i < n; i++)fxym[i] = 6 * f(i-1, i, i+1);for(i = 1; i < n; i++){a[i] = h[i-1] / (h[i] + h[i-1]);c[i] = 1 - a[i];}a[n] = h[n-1] / (h[n-1] + h[n]);cal_m(n);cout<<"\n输出三次样条插值函数:\n";printout(n);cout<<"Do you to have anther try ? y/n :";cin>>ch;}while(ch == 'y' || ch == 'Y');return 0;}void printout(int n){cout<<setprecision(6);for(int i = 0; i < n; i++){cout<<i+1<<": ["<<x[i]<<" , "<<x[i+1]<<"]\n"<<"\t";/*cout<<fxym[i]/(6*h[i])<<" * ("<<x[i+1]<<" - x)^3 + "<<<<" * (x - "<<x[i]<<")^3 + "<<(y[i] - fxym[i]*h[i]*h[i]/6)/h[i]<<" * ("<<x[i+1]<<" - x) + "<<(y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i]<<"(x - "<<x[i]<<")\n";cout<<endl;*/float t = fxym[i]/(6*h[i]);if(t > 0)cout<<t<<"*("<<x[i+1]<<" - x)^3";else cout<<-t<<"*(x - "<<x[i+1]<<")^3";t = fxym[i+1]/(6*h[i]);if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")^3";else cout<<" - "<<-t<<"*(x - "<<x[i]<<")^3";cout<<"\n\t";t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];if(t > 0)cout<<"+ "<<t<<"*("<<x[i+1]<<" - x)";else cout<<"- "<<-t<<"*("<<x[i+1]<<" - x)";t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")";else cout<<" - "<<-t<<"*(x - "<<x[i]<<")";cout<<endl<<endl;}cout<<endl;}。

三次样条插值PPT学习教案

三次样条插值PPT学习教案

x1,x2,…,xn 称为
Sm(x1, x2, …, xn)
第2页/共31页
m=1时,样条函数是分段线性函数; m=2时,是1阶连续可微的分段二次多 项式; 显然, m次样条函数比一般的m次分段插值多 项式的 光滑性 好。 问题: 如何判断一个分段的多项式函数是样条 函数? 设
s(x)
p0 (x), x x1
从而,
a b c 3 3a 2b c 5 6a 2b 8
a 2, b 2, c 3。
第10页/共31页
5.3.2 三 次样条 插值及 其收敛 性(简 介,学 生自学)
有些实际问题中提出的插值问题,要求 插值曲 线具 有较高的光滑性和几何光顺性。 模线员用压铁压住弹性均匀的窄木条 (样条 )的两 端, 强迫样条通过一组已知离散的型值点 。 的形状后,再沿着样条画出所需的曲 线。 形下,该曲线可以由三次样条函数表 示。 插值不仅具有较好的收敛性和稳定性 ,而且 其光滑 性也 较高,因此,样条函数成为了重要的 插值工 具。
进一步可得,
光滑因子
pj (x) pj1(x) qj (x)
第4页/共31页
j 1, 2, , n
对于满足上述性质的如下形式的分段 m次多 项式s(x) ,
p0 (x), x x1
s(x)
p1(x), x1 x x2
pj (x), xj x xj1
pn (x), xn x
pj (x) Pm( j 0,1,...,n)
a 11 a 2 ,
b 1 3 b 2 , c 3。
2)由连续性,应有
ax3 bx2 cx 1 x3 x2
即,
x 1
x 1
a b c 1 13 12 2 a b c 3

三次样条插值

三次样条插值

2012-2013(1)专业课程实践论文三次样条插值樊旭,0818180212,R数学08-2班当插值节点很多时,插值多项式的次数就会很高,这不仅增大了计算量,还会影响结果的精确度。

虽然可以采用原始的插值法,但是主要缺点就是分段接头处不光滑,导数不连续。

因此构想了这样的函数,既能分段的低次插值,又能保证接头处光滑,就产生了三次样条插值。

三次样条插值应该满足的条件:1:满足定义2:S(f x -0)=S(f x +0), S '(f x -0)= S '(f x +0),S ''( f x -0)= S ''(f x +0) 3:边界条件:(1) S '(n x )=0f 'S '(n x )=n f '(2) S ''(0x )=0f 'S ''(n x )=n f '(3) S ''(0x )=S ''(n x )=04:周期函数时:S(0x +0)=S(n x -0),S '(0x +0)=S '(n x -0),S ''(0x +0)=S ''(n x -0),计算方程:1:三转角方程:S(x)=[][][][]11221112212233)()()(2)()(2)(++++++--+--+-+-+-+-j j j j j j j j j j j j j j j j j j m h x x x x m h x x x x y h x x h x x y h x x h x x2:三弯矩方程:S(x)=()j j j j i j j j j i j j J j j j h x x h M y h x x h M y h x x M h x x M -⎪⎭⎫ ⎝⎛-+-⎪⎭⎫ ⎝⎛-+-+-+++++666)(6221113131 以1为例用VC++程序编程就求解NN输入n,xi,yi,和待估点xx ,边界条件。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

摘要本文细致的讲解了三次样条插值函数的产生及在实际中解决的问题,通过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-=+为任意常数。

又根据样条插值函数定义中的条件(3),即k k k y x f x s ==)()(111)()(+++==k k k y x f x s可以确定k a 与k b 为k a =)(611k k k k k k m m hh y y ---++62k kk k h m y b -= (1.04)将式(1.04)中k a 与k b 的值代入表达式(1.03后,就可以得到样条插值函数)(x s 在区间],[1+k k x x 上的表达式为6))((66)(6)()(2113131k kk k k k k k k k k k k k k k h m y x x m m h h y y h x x m h x x m x s -+----+-+-=++++ (1.05)其中k m 与1+k m 分别为区间],[1+k k x x 两端点处的二阶导数值。

由此可以看处,只要能确定各点处的二阶导数值k m )1....,1,0(-=n k ,则子渠道间上的三次样条插值函数)(x s 也确定了。

)(x s 在区间[a,b]上的一阶导数连续,在各节点的左右两子区间上的s(x)虽然不同,但在连接点处的导数存在,即在连接点处的左右导数相等,有)0()0(''+=-k k x s x s(1.06)为了利用条件(2.18),在x 属于],[]1+k k x x 时,县求)('x s 为)(62)(2)()(112121'k k kk k k k k k k k k m m h h y y h x x m h x x m x s ---+-+--=++++(1.07)当x 属于],[1+k k x x 时,11'63)0(++---=+k k k k k k k k m hm h h y y x s(1.08)整理得:)1,.....,2,1(2)1(11-==++-+-n k m m m k k k k k λμμ(1.09)其中1-+=k k kk h h h μ)(61111--+----+=k k k k k k k k k h y y h y y h h λ2 实际问题举例分析例:2.1计算思想设'(0)s m =先求分段Hermite 插值.在[-1,1]上,易知二次函数()22p x x x=+满足(1)1,(0)0,p p -=-= '(1)0p -=为使得三次函数S(x)满足(1)(1)1s p -=-=-,(0)(0)0S p ==,''(1)(1)0,S p -=-=令2()()(1)S x p x c x x =++有'2()22[2(1)(1)]S x x c x x x =+++++.由'(0)2S c m =+=得2c m =-,从而22()2(2)(1)S x x x m x x =++-+在[0,1]上,易得2()(1)p x cx x +-有'2()43[(1)2(1)]S x x c x x x =-++-+-显然(0)(0)0S p ==及(1)(1)1S p ==及''(1)(1)1S p ==-又由'(0)3S c m =+=有3c m =- 从而22()23(3)(1)S x x x m x x =-++-- 求导得2'22(2)[2(1)(1)],10()4(3)[4(1)2],0 1.x m x x x x S x m x x x ⎧++-+++-≤≤=⎨-+--+≤≤⎩最后有''''7(00)(00)24(2)44(3)4S S M m m -=+⇒+-=---⇒=得3232117,10424()517,01424x x x x S x x x x x ⎧-++-≤≤⎪⎪=⎨⎪-++≤≤⎪⎩2.2主要程序代码及命令Function m=naspline(x,y,dy0,dyn,xx)%用途:三样条插值(一阶导数边界条件)%格式:m=naspline(x,y,dy0,xx) x为节点向量,y为数据,dy0,dyn为左右%两短点的一介倒数值,如果xx 缺省,则输出各节点的一阶导数值,m 为xx %的三样条插值n=length(x)-1;h=diff(x);lemda=h(2:n)./(h(1:n-1)+h(2:n));mu=1-lemda;g=3*(lemda.*diff(y(1:n))./h(1:n-1)+mu.*diff(y(2:n+1))./h(2:n));g(1)=g(1)-lemda(1)*dy0;g(n-1)= g(n-1)-mu(n-1)*dyn;%求解三对角方程dy=nachase(lemda,2*ones(1:n-1),mu,g);m=[dy0;dy;dyn];if nargin>=5s=zeros(size(xx));for i=1:nif i==1,kk=find(xx<=x(2));elseif i==nkk=find(xx>x(n));elsekk=find(xx>x(i)&xx<=x(i+1));endxbar=(xx(kk)-x(i))/h(i);s(kk=alpha0(xbar)*y(i)+alpha1(xbar)*y(i+1)+…h(i)*beta0(xbar)*m(i)+h(i)*beta1(xbar)*m(i+1);endm=s;endfunction x=nachase(a,b,c,d)n=length(a);for k=2:nb(k)=b(k)-a(k)/b(k-1)*c(k-1);d(k)=d(k)-a(k)/b(k-1)*d(k-1);endx(n)=d(n)/b(n);for k=n-1:-1:1x(k)=(d(k)-c(k)*x(k+1))/b(k);endx=x( : );function y=alpha0(x)y=2*x.^3-3*x.^2+1;function y=alpha1(x)y=-2*x.^3+3*x.^2;function y=beta0(x)y=x.^3-2*x.^+x;function y=beta1(x)y=x.^3-x.^2;2.3结果演示naspline([-1 0 1],[-1 0 1],0,-1) %输入m的值ans=1.7500-1.0000naspline([-1 0 1],[-1 0 1],0,-1,-1:0.25:1)ans=-1.0000 -0.9258 -0.7188 -0.4023 0 -2.2444 -1.631 -1.10981 -0.7500总结随着计算机技术及硬件设施的不断发展,计算机语言的演化从最开始的机器语言到汇编语言,最后到支持面向对象技术的面向对象语言。

这就要求提供的计算方法也不断发展。

同时在实践也给旧的插值逼近方法不断提出问题,这就要求新的插值逼近方法要更复杂,误差精度更高,同时解决更多方面的问题。

插值逼近方法的发展也需要新的理论指导。

自然辩证法的科学理论中提到科学“范式”概括了插值逼近法得法的发展过程.辩证唯物主义自然观、自然科学发展过程及其规律,分析与综合、归纳法与演绎法、想象和类比等科学逻辑思维方法的应用都在插值逼近理论的发展过程中起到了重要作用。

用科学的逻辑思维方法认识事物才会清楚的了解其过去、现在和未来.计算数学中的插值逼近方法发展同样遵循着科学技术、科学理论发展的一般规律.以自然辩证法的观点来分析.有助于我们更加深入地认识插值逼近以及整个数值计算方法发展的历史、现状和趋势。

插值逼近方法及数学理论的进一步发展也必将为自然辩证法的发展提供基础。

本次的课程设计的整个过程让我认识到了基础知识的欠缺,我通过查阅大量的资料,重新阅读关于计算方法、数学分析、高等代数、MATLAB方面的书籍,从根源上了解,三样条插值的由来和应用范围等,是我受益良多。

切切实实的认识的努力学习的重要性。

参考资料[1]陈辉,李文宇,张传芳数值计算方法哈尔滨哈尔滨工业大学出版社,2009[2]李庆扬,易大义,王能超. 现代数值分析. 北京高等教育出版社,1995[3]刘春风,何亚丽,应用数值分析北京冶金工业出版社 2005[4]郝红伟,MATLAB 6,北京,中国电力出版社,2001[5]姜健飞,胡良剑,数值分析及其MATLAB实验,科学出版社,2004[6]薛毅,数值分析实验,北京工业大学出版社,2005。

相关文档
最新文档