三次样条函数

计算方法实验报告

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

h[k] = x[k+1] - x[k];

/*构建Am=f中的行列式A*/

for(k=1; k

{

rr[k] = h[k]/(h[k]+h[k-1]); /*计算λ[k]*/

uu[k] = h[k-1]/(h[k]+h[k-1]); /*计算μ[k]*/

gg[k] = 3 * (uu[k]*(y[k+1]-y[k])/h[k] + rr[k]*(y[k]-y[k-1])/h[k-1]); /*计算g[k]*/

}

/* g[1] = λ[1]*f＇[0] */

/* g[n-1] = μ[n-1]*f＇[n] */

gg[1] = gg[1] - rr[1]/x[0];

gg[n-1] = gg[n-1] -uu[n-1]/x[n];

/*作A=LU分解*/

aa[1] = 2;

bb[1] = uu[1]/aa[1];

for(k=2; k

{

aa[k] = (float)2.0 - rr[k]*bb[k-1];

bb[k] = uu[k]/aa[k];

}

/*通过LT＝G计算中间矩阵T*/

tt[1] = gg[1]/aa[1];

for(k=2; k

tt[k] = (gg[k]-rr[k]*tt[k-1])/aa[k];

/*根据第一类边界条件所得*/

mm[n] = 1/x[n];

mm[0] = 1/x[0];

/*通过UM=T计算M*/

mm[n-1] = tt[n-1];

for(k=n-2; k>=1; k--)

mm[k] = tt[k]-bb[k]*mm[k+1];

do{

printf("输入插值点：");

scanf("%f",&u);

/*判断u属于哪个区间x[k],x[k+1]确定k*/

for(k=0; k

{

if( u>=x[k] && u<=x[k+1] )

{

kv = k;

break;

}

}

/*代入u计算值*/

s = (1-2*(u-x[kv])/(-h[kv])) * ((u-x[kv+1])/(-h[kv])) * ((u-x[kv+1])/(-h[kv])) * y[kv];

s = s + (1-2*(u-x[kv+1])/h[kv]) * ((u-x[kv])/h[kv]) * ((u-x[kv])/h[kv]) * y[kv+1];

s = s + (u-x[kv])*((u-x[kv+1])/(-h[kv]))*((u-x[kv+1])/(-h[kv])) * mm[kv];

s = s + (u-x[kv+1])*((u-x[kv])/(-h[kv]))*((u-x[kv])/(-h[kv])) * mm[kv+1];

printf("与%f对应的结点值为%f\n\n",u,s);

printf("还要继续插值吗?(1/0): ");

scanf("%d",&chose);

}while( chose != 0 );

}

void main()