ex4.cpp

方程的数值计算方法

//第四章，3次样条插值
#include<iostream.h>
#include<stdlib.h>
#define N 20
double result[N];//存放M矩阵

double ChaShang(double x[N],double fx[N],double j)//j=1,M-1(n-2)
{
	double f1;
	double f2;
	f1=(fx[j+1]-fx[j])/(x[j+1]-x[j]);
	f2=(fx[j]-fx[j-1])/(x[j]-x[j-1]);
	if(x[j+1]!=x[j-1])
		return (f1-f2)/(x[j+1]-x[j-1]);
	else
	{
		cout<<"计算插商时出现错误！";
		exit(0);
	}
}
void CalH(double x[N],double h[N],int n)//计算两个子变量之间的距离
{
	for(int i=0;i<n-1;i++)
		h[i]=x[i+1]-x[i];
}
void Calpara(double h[N],double para1[N],double para2[N],int n)//计算解方程组中用到的参数
{
	for(int i=1;i<n-1;i++)
	{
		para1[i]=h[i-1]/(h[i-1]+h[i]);
		para2[i]=1-para1[i];
	}
}
void SetValue(double matrix[N][N+1],double para1[N],double para2[N],double b[N],int n)
{
	matrix[0][n]=b[0];
	for(int i=1;i<n-1;i++)
	{
		matrix[i][i-1]	=para1[i];
		matrix[i][i]	=2;
		matrix[i][i+1]	=para2[i];
		matrix[i][n]	=b[i];
	}
	matrix[n-1][n-2]	=1;
	matrix[n-1][n-1]	=2;
	matrix[n-1][n]		=b[n-1];
}
void CalB(double b[N],double x[N],double fx[N],double h[N],int n)
{
	for(int j=1;j<n-2;j++)
	{
		b[j]=6*(ChaShang(x[j],fx[j],x[j+1],fx[j+1])-ChaShang(x[j-1],fx[j-1],x[j],fx[j]))/(h[j]-h[j-1]);
	}
}
//消元运算
void elimination(double a[N][N+1],int n)	
{
	int i=0,j=0;
	int row=0;
	double value=a[0][0];

	for(i;i<n;i++)							//大循环，实现所有行的消元
	{
		for(j=i+1;j<n;j++)		//row
		{
			value=a[j][i]/a[i][i];
			for(int k=i;k<=n;k++)//每行的元素
			{
				a[j][k]+=a[i][k]*value*(-1);
			}
		}
	}
}

//回代运算
void ReCalculation(double a[N][N+1],int n)
{
	double value=0;
	if(a[n-1][n-1]==0)			//首先确定分母不为零
	{
		cout<<"方程组无解";
		return;
	}
	result[n-1]=a[n-1][n]/a[n-1][n-1];
	
	for(int i=n-2;i>=0;i--)
	{
		value=0;
		for(int j=i+1;j<n;j++)
		{
			value+=a[i][j]*result[j];
		}
		result[i]=(a[i][n]-value)/a[i][i];
	}
}
void main()
{
	double x[N];
	double fx[N];
	double h[N],para1[N]={0,0};
	double para2[N]={0,0};
	double b[N];
	int n;
	double matrix[N][N+1]={{2,1}};
	double boundary1,boundary2;

	cout<<"请输入节点个数：";
	cin>>n;			//如果不是一个整数，则用舍去法取整 n=M+1,0，1...M是一维数组最大的下标

	cout<<"\n请输入节点自变量的值：\n";
	for(int i=0;i<n;i++)
	{
		cin>>x[i];
	}

	cout<<"请输入"<<n<<"个节点的函数值：\n";
	for(i=0;i<n;i++)
	{
		cin>>fx[i];
	}
	cout<<"请输入边界条件，Y("<<x[0]<<")'=";
	cin>>boundary1;
	cout<<"Y("<<x[n-1]<<")'=";
	cin>>boundary2;
	CalH(x,h,n);
	b[0]	=6*(fx[1]-boundary1)/h[0];
	b[n-1]	=6*(boundary2-fx[n-1])/h[n-1];
	CalH(x,h,n);
	Calpara(h,para1,para1,n);
	CalB(b,x,fx,h,n);
	SetValue(matrix,para1,para2,b,n);
	elimination(matrix,n);
	ReCalculation(matrix,n);
	for(i=0;i<n;i++)
	{
		cout<<result[i];
	}

}