样条插值.cpp

这是一些c++例程

// 样条插值.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "stdio.h"
#include "math.h"

#define SAMS 19
#define WIDTH 78
#define HEIGHT 22


class CSample
{
public:
	char screen[HEIGHT][WIDTH];			//成员变量，表示文本图像
	double x[SAMS], y[SAMS];			//成员变量，分别表示样点的横坐标与纵坐标
	double m[SAMS];					//成员变量，表示样点处的导数值。

	CSample();						//构造函数，初始化文本图像
	void Tribble();					//成员函数，进行三次样条插值。
	void Add(int, double);				//成员函数，用来向screen数组中加入星号
	void Show();						//成员函数，用来显示文本图
};

CSample::CSample()                                //构造函数，初始化文本图像
{
	for(int i=0;i<HEIGHT;i++)
		for(int j=0;j<WIDTH;j++)
			screen[i][j]=' ';				//开始时全部置为空格
}

void CSample::Show()						//成员函数，用来显示文本图
{
	for(int i=0;i<HEIGHT;i++)
	{
		for(int j=0;j<WIDTH;j++)
			printf("%c",screen[i][j]);			//逐字输出文本图
		printf("\n");
	}
}

void CSample::Add(int a,double b)				//成员函数，用来向screen数组中加入星号
{
	int x0=0,y0=21;						//机翼曲线的原点位置
	screen[y0-int(b/1.7)][x0+a/7]='*';				//把机翼曲线图像对应的点置为星号
}

void CSample::Tribble()
{
	int n=SAMS-1;
	
	for(int _x=1;_x<521;_x++)						//_x表示横坐标
	{
		double u=_x;							//将横坐标转为double型参与运算
		double res=0;							//res代表纵坐标插值运算结果
		double *lanmda=new double[n];
		double *miu=new double[n];
		double *q=new double[n];
		for(int i=1;i<n;i++)
		{
			lanmda[i]=(x[i+1]-x[i])/(x[i+1]-x[i-1]);
			miu[i]=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
			q[i]=3*(lanmda[i]*(y[i]-y[i-1])/(x[i]-x[i-1])+miu[i]*(y[i+1]-y[i])/(x[i+1]-x[i]));
		}
		q[1]-=lanmda[1]*m[0];
		q[n-1]-=miu[n-1]*m[n];
		
		double *r=new double[n];
		double *p=new double[n];
		double *g=new double[n];
		p[1]=2; g[1]=q[1]/p[1];
		for(i=2;i<n;i++)
		{
			r[i-1]=miu[i-1]/p[i-1];
			p[i]=2-lanmda[i]*r[i-1];
			g[i]=(q[i]-lanmda[i]*g[i-1])/p[i];
		}
		m[n-1]=g[n-1];
		for(i=n-2;i>=1;i--)
			m[i]=g[i]-r[i]*m[i+1];
		delete[] lanmda;
		delete[] miu;
		delete[] q;
		delete[] r;
		delete[] p;
		delete[] g;
		
		int j;
		double tmp1,tmp2;
		for(i=0;i<=n;i++)
			if (x[i]<=u && u<=x[i+1]) j=i;
			if (j<n) 
			{
				tmp1=(u-x[j+1])/(x[j]-x[j+1]);
				tmp2=(u-x[j])/(x[j+1]-x[j]);
				
				res+=y[j]*tmp1*tmp1*(1+2*(u-x[j])/(x[j+1]-x[j]));
				res+=m[j]*tmp1*tmp1*(u-x[j]);
				j++;
				res+=y[j]*tmp2*tmp2*(1+2*(u-x[j])/(x[j-1]-x[j]));
				res+=m[j]*tmp2*tmp2*(u-x[j]);
			}
			if (_x%7==0) Add(_x,res);				//返回方式：向screen的对应为置加入星号
	}
}

void main()
{
	CSample test;							//实例化CSample类
	//给出各样点坐标
	test.x[0]=0.520;		test.y[0]=5.288;
	test.x[1]=3.1;			test.y[1]=9.4;
	test.x[2]=8.0;			test.y[2]=13.84;
	test.x[3]=17.95;		test.y[3]=20.20;
	test.x[4]=28.65;		test.y[4]=24.90;
	test.x[5]=39.62;		test.y[5]=28.44;
	test.x[6]=50.65;		test.y[6]=31.10;
	test.x[7]=78.0;			test.y[7]=35.0;
	test.x[8]=104.6;		test.y[8]=36.9;
	test.x[9]=156.6;		test.y[9]=36.6;
	test.x[10]=208.6;		test.y[10]=34.6;
	test.x[11]=260.7;		test.y[11]=31.0;
	test.x[12]=312.50;		test.y[12]=26.34;
	test.x[13]=364.4;		test.y[13]=20.9;
	test.x[14]=416.3;		test.y[14]=14.8;
	test.x[15]=468.0;		test.y[15]=7.8;
	test.x[16]=494.0;		test.y[16]=3.7;
	test.x[17]=507.0;		test.y[17]=1.5;
	test.x[18]=520.0;		test.y[18]=0.2;

	//给出两端点的导数
	test.m[0]=1.86548; test.m[SAMS-1]=-0.046115;

	test.Tribble();								//进行运算
	test.Show();								//输出结果
}