lathinsertvalue.cpp

常用数值算法的C++源码

#include "stdio.h"
#include "math.h"
                        
double *h;                        //步长序列
double *dif2;                     //2阶差商序列
double *d,*mu,*lamed,*a,*m;       //第一型3次样条插值函数矩阵方程元素       

//高斯消去法
void Gauss(double *a,int n)
{
	int    col,N=n+1;
	double max,temp,sum;
	for(int k=0;k<n-1;k++)
	{
		max=fabs(a[k*N+k]);
		col=k;
		for(int i=k+1;i<n;i++)
		{
			temp=fabs(a[i*N+k]);
			if(temp>max)
			{
				max=temp;
				col=i;
			}
		}
		if(col!=k)
		{
			for(i=k;i<N;i++)
			{
				temp=a[col*N+i];
				a[col*N+i]=a[k*N+i];
				a[k*N+i]=temp;
			}
		}
		for(i=k+1;i<n;i++)
		{
			if(a[i*N+k]!=0)
			{
				for(int j=k;j<N;j++)
					m[j]=a[i*N+j]-a[k*N+j]*a[i*N+k]/a[k*N+k];
				for(j=k;j<N;j++)
					a[i*N+j]=m[j];
			}
		}
	}
    m[n-1]=a[(n-1)*N+n]/a[(n-1)*N+(n-1)];
	for(int i=n-2;i>=0;i--)
	{
		sum=0;
		for(int j=i+1;j<n;j++)
			sum=sum+a[i*N+j]*m[j];
		m[i]=(a[i*N+n]-sum)/a[i*N+i];
	}
}

//第一型3次样条插值函数
void LathInsertValue(double *x,double *y,double dev_0,double dev_n,int n)
{
	//计算步长序列和2阶差商序列
	for(int i=0;i<n;i++)
	{
		h[i]    = x[i+1]-x[i];
		dif2[i] = (y[i+1]-y[i])/h[i];
	}
	//计算系数矩阵元素
	for(i=1;i<n;i++)
	{
		mu[i]    = h[i-1]/(h[i-1]+h[i]);
		lamed[i] = 1-mu[i];
		d[i]     = 6*(dif2[i]-dif2[i-1])/(h[i-1]+h[i]);
	}
	d[0] = 6*(dif2[0]-dev_0)/h[0];
	d[n] = 6*(dev_n-dif2[n-1])/h[n-1];
	//生成系数矩阵
	for(i=1;i<n;i++)
	{
		a[i*(n+2)+(i-1)] = mu[i];
		a[i*(n+2)+i]     = 2;
		a[i*(n+2)+(i+1)] = lamed[i];
		a[i*(n+2)+n+1]   = d[i];
	}
	a[0]             = 2;
	a[1]             = 1;
	a[n+1]           = d[0];
	a[n*(n+2)+(n-1)] = 1;
	a[n*(n+2)+n]     = 2;
    a[n*(n+2)+n+1]   = d[n];
	//追赶法解矩阵方程
	Gauss(a,n+1);
}

void main()
{
	int    n     = 10;
	double dev_0 = 0.8;
	double dev_n = 0.2;
	double x[11] = {0,1,2,3,4,5,6,7,8,9,10};
	double y[11] = {2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80};
	double h_g[10],dif2_g[10],mu_g[10],lamed_g[10],d_g[11],m_g[11],a_g[11*12],s[11];
	FILE *fp=fopen("LathInsertValue.txt","w");
	h       = h_g;
	dif2    = dif2_g;
	mu      = mu_g;
	lamed   = lamed_g;
	d       = d_g;
	m       = m_g;
	a       = a_g;
	for(int i=0;i<n+1;i++)
		for(int j=0;j<n+1;j++) a[i*(n+2)+j] = 0;
	LathInsertValue(x,y,dev_0,dev_n,n);
	for(i=0;i<n;i++)
	{
		s[i] = y[i]+(dif2[i]-(m[i]/3+m[i+1]/6)*h[i])*0.5+0.5*m[i]*pow(0.5,2)+(m[i+1]-m[i])*pow(0.5,3)/h[i]/6;
		fprintf(fp,"S(%.1f) = %f\n",i+0.5,s[i]);
		printf("S(%.1f) = %f\n",i+0.5,s[i]);
	}
	printf("\n运行结果保存至文件LathInsertValue.txt\n");
}