样条插值.cpp

分析研究基于三种边界条件的样条插值算法并用程序实现其中一种边界条件

/*
注:所用的数据表如下:
x 0 70 130 210 337 578 776 1012 1142 1462 1841
y 0 57 78   103 135 182 214 244   256   272    275
边界条件为:
(1)y'(0)=1,y'(1841)=0;
(2)y"(0)=0,y"(1841)=0.
求解在以上边界条件下在x[k]=50k(k=1,2...,36)处的函数值

程序源代码*/

/*利用三次样条插值法求曲线在某个点处的函数值(本题用第二个边界条件)*/
#include<stdio.h>
#include<math.h>
#define N 11            //N表示节点的个数

void main()
{
	double x[N]={0,70,130,210,337,578,776,1012,1142,1462,1841};
	double y[N]={0,57,78,103,135,182,214,244,256,272,275};
	double h[N],a[N],b[N],A[N],B[N],m[N],s[37],xx[37]; //s[N]表示曲线,其中数组xx[N]表示自变量x[k]
	int i,k;   
	h[0] = x[1]-x[0];       //初始化h0,a0,b0,A0,B0
	a[0] = 1;
	b[0] = 3*(y[1]-y[0])/h[0];
	A[0] = -a[0]/2;
	B[0] = b[0]/2;
	for(i=0;i<N;i++)       //求hi
		h[i]=x[i+1]-x[i];
	
	for(i=1;i<N-1;i++)//求ai,bi
	{   
		a[i]=h[i-1]/(h[i-1]+h[i]);
		b[i] = 3 * ((1-a[i])*(y[i]-y[i-1])/h[i-1] + a[i]*(y[i+1]-y[i])/h[i]);
	}
	
	for(i=1;i<N-1;i++) 
	{ //求Ai,Bi
		A[i] = -a[i]/(2+(1-a[i])*A[i-1]);
		B[i] = (b[i]-(1-a[i])*B[i-1])/(2+(1-a[i])*A[i-1]);
	}
	
	m[N-1]=(b[N-1]-(1-a[N-1])*B[N-2])/(2+(1-a[N-1])*A[N-2]);//求Mn的值

	
	for(i=N-2;i>=0;i--)        //求m0,m1,-----mn-1的值
		m[i] = A[i] * m[i+1] + B[i];
	
	for(k=1;k<=36;k++) 
	{       //求曲线在x[k]处的函数值
		xx[k] = 50 * k;
		
		for(i=0;i<N;i++)      //找出x[k]所在的区间
			if(xx[k] >= x[i] && xx[k] <= x[i+1])
				//s[k]即为x[k]所在区间的三次样条插值函数,以下即为求在x[k]处的函数值
		s[k] = (1+2*(xx[k]-x[i])/(x[i+1]-x[i]))*pow((xx[k]-x[i+1])/(x[i]-x[i+1]),2)*y[i] +
       (1+2*(xx[k]-x[i+1])/(x[i]-x[i+1]))*pow((xx[k]-x[i])/(x[i+1]-x[i]),2)*y[i+1] +
       (xx[k]-x[i])*pow((xx[k]-x[i+1])/(x[i]-x[i+1]),2)*m[i] +
       (xx[k]-x[i+1])*pow((xx[k]-x[i])/(x[i+1]-x[i]),2)*m[i+1] ;
	}
	
	printf("所求的外形曲线在x[k]=50*k(k=1,2,...,36)处的函数值分别为:\n");
	for(k=1;k<=36;k++)        //输出在x[k]处的函数值
		printf("s(%4d) = %13.8f\n",50*k,s[k]);
	printf("\n");
	
	printf("m[i]的值分别为:\n");
	for(i=0;i<N;i++)
		printf("m(%2d) = %8f\n",i,m[i]);
	printf("\n");
}