数值计算之插值法

前段时间，接触了一下数值计算，课程下来，数学水平没有太多见长，到时复习了一下C++和一些数据结构的知识。自己动手敲了下代码，实现了几个常见的数值公式。插值法，现在想来，其实是在通过有限的点或者说是采集的样本点，来拟合成可能的函数关系式。下面是对一组数据的拟合的例子。

#include "stdafx.h"
#include<iostream>
#include<cmath>
#include<math.h>
using namespace std;
const double e=2.718281828459045;
double getSum(int m,int f1,int f2);//计算法方程等式左边矩阵中各个值
double getYSum(int n,int f);//计算法方程等式右边矩阵中的各个值
double getSquareError(int n);//计算平方误差公式的第一项的值
const int M=50;
double X[M];
double Y[M];
double W[M];
double Func[M][M];//存放各个拟合函数
double xMatrix[M][M];//存放法方程等式左边的矩阵
double yMatrix[M];//存放法方程等式右边的矩阵

class Matrix{
	int size;//矩阵的大小
	const static int N=100;
	double matrix[N][N];//增广矩阵
	double solution[N];//方程的解
public:
	Matrix(int k);//构造系数矩阵
	void disPlay();//显示增广矩阵
	void setAugmetMatrix(double *BMatrix);//生成增广矩阵
	bool solveMatrix();//解矩阵，其中有消元和回代的过程
	int getSize();//获取矩阵的大小
	double* getSolution();//获得解的数组
	double* getFirstMatrix();//获取矩阵第一行第一列元素
	int getMaxMainElement(int currentCol);//获取最大主元所在的行
	void exchangTwoRow(int row1,int row2);//交互矩阵中的两行
};
bool Matrix::solveMatrix(){
#pragma region 矩阵消元
	for (int i = 1; i <= size; i++)//最外层循环遍历整个矩阵
	{
		if (getMaxMainElement(i)>0 &&i!=getMaxMainElement(i))
		{
			exchangTwoRow(getMaxMainElement(i),i);//交换行
		}
		cout<<"第"<<i<<"变换"<<endl;
		disPlay();
		cout<<endl;
		for (int j = i+1; j <= size; j++)//第二层循环遍历与上一层相邻的层
		{
			double temp=matrix[j][i];//将在第k次消去的主元的第一个元素赋给temp临时变量
			matrix[j][i]=0;//将第一元素置0
			for (int k = i+1; k <=size+1; k++)//第三层循环遍历矩阵的列
			{
				matrix[j][k]=matrix[j][k]-(temp/matrix[i][i])*matrix[i][k];//根据矩阵消去公式依次肖元，成为上三角矩阵
			}
		}

	}  
#pragma endregion
#pragma region 回代
	if (matrix[size][size]!=0)//预先判断消去后的矩阵是否有确定解
	{
		double temp;
		solution[1]=matrix[size][size+1]/matrix[size][size];//解得xn存入solution数组
		for (int i =2; i <=size; i++)
		{
			double A=matrix[size-i+1][size+1];//将bi赋给变量A
			for (int k = 1,j=1; k <i; j++,k++)
			{
				temp=A-solution[j]*matrix[size-i+1][size-k+1];//做减法
				A=temp;//赋值使上面等式可连减
			}
			solution[i]=temp/matrix[size-i+1][size-i+1];//除以主对角线上的系数，得到xi
		}
		return true;
	}  
#pragma endregion

	return false;
}
Matrix::Matrix(int n){
	size=n;
	cout<<"请输入方程的系数矩阵"<<endl;
	for(int i=1;i<=n;i++)
		for (int j = 1; j <=n; j++)
		{
			matrix[i][j]=xMatrix[i-1][j-1];//录入矩阵
		}
}
double* Matrix::getFirstMatrix(){
	double*p=& matrix[1][1];
	return p;
}
double* Matrix::getSolution(){
	return solution;
}
int Matrix::getSize(){
	return size;
}
void Matrix::disPlay(){
	for (int k = 1; k <= size; k++)
	{
		for (int m = 1; m <=size+1; m++)
		{
			if (m==size+1)
			{
				cout<<matrix[k][m]<<endl;//
			}
			else
			{
				cout<<matrix[k][m]<<"\t";
			}
		}
	}
}
void Matrix::setAugmetMatrix(double*BMatrix){
	for (int i = 1; i <=size; i++)
	{
		matrix[i][size+1]=BMatrix[i];//
	}
}
int Matrix::getMaxMainElement(int currentCol){
	double temp=0;
	int maxRowMain=0;
	for (int i = currentCol; i <= size; i++)
	{
		if (temp<matrix[i][currentCol])
		{
			temp=matrix[i][currentCol];
			maxRowMain=i;
		}
	}
	return maxRowMain;
}
void Matrix::exchangTwoRow(int row1,int row2){
	double tempMatrix=0;
	for (int i = 1; i <=size+1; i++)
	{
		tempMatrix=matrix[row1][i];
		matrix[row1][i]=matrix[row2][i];
		matrix[row2][i]=tempMatrix;
	}
}
void _tmain(int argc, _TCHAR* argv[])
{
	int num;
	int n;
	cout<<"请输入点的个数：";
	cin>>n;
	cout<<"输入拟合函数个数:";
	cin>>num;
	//输入点的值
	cout<<"输入各个点的值"<<endl;
	for (int i = 0; i < n; i++)
	{
		cout<<"x"<<i<<"=";
		double x;
		cin>>x;
		X[i]=x;

		cout<<"y"<<i<<"=";
		double y;
		cin>>y;
		Y[i]=y;

		cout<<"w"<<i<<"=";
		double w;
		cin>>w;
		W[i]=w;
		cout<<endl;
	}
	//S(x)=a*ln(x)+b*cos(x)+c*e^x
	for (int j = 0; j < n; j++)
	{
		Func[0][j]=log(X[j]);
		Func[1][j]=cos(X[j]);
		Func[2][j]=pow(e,X[j]);
	}


	for (int i = 0; i < num; i++)
	{
		for (int j = i; j < num; j++)
		{
			xMatrix[i][j]+=getSum(n,i,j);
			if (i!=j)
			{
				xMatrix[j][i]=xMatrix[i][j];
			}
		}
	}

	for (int i = 1; i <= num; i++)
	{
		yMatrix[i]=getYSum(n,i);
	}
		Matrix m(num);
	m.setAugmetMatrix(yMatrix);
	if(m.solveMatrix()){
		m.disPlay();//该方法有测试矩阵运算是否正确之用
		double* solute=m.getSolution();
		cout<<"方程的解为："<<endl;
		for (int i =0; i < m.getSize(); i++)
		{
			char a='a'+i;
			cout<<a<<"="<<solute[m.getSize()-i]<<endl;
		}
		//计算平方误差
		double sError=getSquareError(n)-yMatrix[1]*(solute[m.getSize()])-yMatrix[2]*(solute[m.getSize()-1])-yMatrix[3]*(solute[m.getSize()-2]);
		cout<<"平方误差："<<sError<<endl;
	}
	else
	{
		cout<<"此方程无解或无准确解。"<<endl;
	}


}
double getSum(int n,int f1,int f2)
{
	double result=0;
	for (int i = 0; i < n; i++)
	{
		result+=W[i]*Func[f1][i]*Func[f2][i];
	}
	return result;
}

double getYSum(int n,int f)
{
	double result=0;
	for (int i = 0; i < n; i++)
	{
		result+=W[i]*Func[f-1][i]*Y[i];
	}
	return result;
}

double getSquareError(int n)
{
	double result=0;
	for (int i = 0; i < n; i++)
	{
		result+=W[i]*Y[i]*Y[i];
	}
	return result;
}