前段时间,接触了一下数值计算,课程下来,数学水平没有太多见长,到时复习了一下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;
}