mymath.h

数值计算方法中常用方法的程序实现

#ifndef MYMATH_H
#define MYMATH_H
//-------------------------
#include<iostream>
#include<vector>
#include<iomanip>
#include<fstream>
#include<sstream>
#include<cmath>
#include<cassert>
#include<utility>
//--------------------------
const double EPS = 1e-10;


//--------------------------
//3维点坐标类，支持2维点和基本点操作
class Point
{
	double x;
	double y;
	double z;
public:
	Point(double xx=0, double yy=0, double zz=0): x(xx), y(yy), z(zz) {}
	Point(const Point& p) { x=p.x; y=p.y; z=p.z; }
	void setX(double a) { x=a; }
	double getX() const { return x; }
	void setY(double b) { y=b; }
	double getY() const { return y; }
	void setZ(double c) { z=c; }
	double getZ() const { return z; }
	void MoveTo(double a=0, double b=0, double c=0) { x=a; y=b; z=c; }
	Point operator+(const Point& p) const { return Point(x+p.x, y+p.y, z+p.z); } 
	Point& operator+=(const Point& p) { x+=p.x; y+=p.y; z+=p.z; return *this; }
	friend inline std::ostream& operator<<(std::ostream& o, const Point& p)
	{
		return o<<'('<<p.x<<','<<p.y<<','<<p.z<<')'<<'\n';
	}
};
//-----------------------------------
//二次贝赛尔曲线函数
//2次曲线需要3点，首尾2点，因为型线近似看成90度，由首尾2个点的坐标可以得到第三点坐标
void Bezier(const Point& a, const Point& b, std::size_t n, std::vector<Point>& vec)
{
	for(std::size_t i=0; i!=n; ++i)
	{
		double u=i*1.0/(n-1);
		double x=a.getX()*pow(1-u,2.0)+2*u*(1-u)*b.getX()+b.getX()*pow(u,2.0);
		double y=a.getY()*pow(1-u,2.0)+2*u*(1-u)*a.getY()+b.getY()*pow(u,2.0);
		Point p(x,y);
		vec.push_back(p);
	}
}
//--------------------------------------
//一般高斯消去法
//消去过程
void Elim(std::vector<std::vector<double> >& v)
{
	for(std::size_t i=0; i!=v.size()-1; ++i)
		for(std::size_t j=i+1; j!=v.size(); ++j)
		{
			double c=v[j][i]/v[i][i];
			for(std::size_t k=i; k!=v[i].size(); ++k)
				v[j][k]-=c*v[i][k];
		}
}
//回代过程
void Back(std::vector<std::vector<double> >& a, std::vector<double>& b)
{
	std::size_t m=a.size();    
	std::size_t n=a[0].size(); 
	std::size_t t=b.size();    
	b[t-1]=a[m-1][n-1]/a[m-1][m-1];  
	for(int i=m-2; i>=0; --i)          //i必须为int型，不能为unsigned，unsigned 减到-1就绕回了
	{
		double c=a[i][n-1];            
		for(unsigned j=i+1; j!=m; ++j)
			c-=a[i][j]*b[j];                 
		b[i]=c/a[i][i];                  
	}
}
//-------------------------------------
//列主元高斯消去法
int Gauss(std::vector<std::vector<double> >& a, std::vector<double>& b, std::vector<double>& x)
	//系数矩阵为向量a，向量x用来存储计算结果
{
	std::size_t n=a.size();
	std::cout<<std::fixed;
	std::cout<<"高斯列主消元法解线性方程组:\n\n方程组增广矩阵：\n";
	for(std::size_t i=1; i<=n; i++)
	{	
		for(std::size_t j=1; j<=n; j++)
			std::cout<<a[i-1][j-1]<<'\t';
		std::cout<<b[i-1]<<std::endl;
	}

	//循环主体
	for(std::size_t k=1; k<=n-1; k++) 
	{	//首先按列选主元
		double zhu=fabs(a[k-1][k-1]);
		std::size_t zhu_i=k;
		for(std::size_t i=k; i<=n; i++)//i从第k行搜索到第n行
			if(fabs(a[i-1][k-1])>zhu)
			{	
				zhu=fabs(a[i-1][k-1]);
				zhu_i=i;
			}
		if(zhu_i!=k)
		{	//交换第zhu_i行和第k行
			double m=b[k-1];
			b[k-1]=b[zhu_i-1];
			b[zhu_i-1]=m;
			for(std::size_t j=k; j<=n; j++)
			{	
				m=a[k-1][j-1];
				a[k-1][j-1]=a[zhu_i-1][j-1];
				a[zhu_i-1][j-1]=m;
			}
		}		
		if(fabs(a[k-1][k-1])<EPS) 
		{
			std::cout<<"No Solution!"<<std::endl;
			return 0;
		}
		for(std::size_t i=k+1; i<=n; i++)
		{	
			double m=-a[i-1][k-1]/a[k-1][k-1];
			for(std::size_t j=k;j<=n;j++)
				a[i-1][j-1]+=m*a[k-1][j-1];
			b[i-1]+=m*b[k-1];
		}
	}
	//循环主体结束
	if(fabs(a[n-1][n-1])<EPS) 
	{
		std::cout<<"No Solution!\n";
		return 0;
	}
	//消元后的上三角增广矩阵
	std::cout<<"\n消元后方程组上三角增广矩阵：\n";
	for(std::size_t i=1; i<=n; i++)
	{	
		for(std::size_t j=1; j<=n; j++)
			std::cout<<a[i-1][j-1]<<'\t';
		std::cout<<b[i-1]<<std::endl;
	}	
	//开始上三角矩阵 回代过程  解存储在x中
	x[n-1]=b[n-1]/a[n-1][n-1];
	for(int i=n-1; i>=1; i--)
	{	
		double m=0;
		for(std::size_t j=i+1;j<=n;j++)
			m+=a[i-1][j-1]*x[j-1];
		x[i-1]=(b[i-1]-m)/a[i-1][i-1];
	}
	//打印解
	std::cout<<"\n解：\n";
	for(std::size_t i=1;i<=n;i++)
		std::cout<<'x'<<i<<": "<<x[i-1]<<std::endl;
	return 1;
}
//-------------------------------------
//追赶法解三对角方程组 (Ax=b)
//从左至右，a对应3对角线第一斜行元素，b为对角线元素，c为第三斜行元素，abc构成了系数矩阵A
//d为方程Ax=b中对应的b，x为所求元素
void Chase(std::vector<double>& a, std::vector<double>& b, std::vector<double>& c, 
		  std::vector<double>& d, std::vector<double>& x)
{
	std::size_t n = b.size();
	std::vector<double> l(n,0);
	std::vector<double> r(n-1,0);
	std::vector<double> y(n,0);
	//克罗特分解过程（追）
	l[0]=b[0];
	y[0]=d[0]/l[0];
	for(std::size_t i=1; i!=n; ++i)
	{
		r[i-1]=c[i-1]/l[i-1];
		l[i]=b[i]-a[i-1]*r[i-1];
		y[i]=(d[i]-a[i-1]*y[i-1])/l[i];
	}
	//回代过程（赶）
	x[n-1]=y[n-1];
	for(int i=n-2; i>=0; --i)
	{
		x[i]=y[i]-r[i]*x[i+1];
	}
}
//-------------------------------------
//杜里特尔分解方程组 (Ax=b)
//a为系数矩阵，b为方程中的b，x为所求
//l为下三角阵，r为上三角阵，y是Rx=y中的y
void Dltr(std::vector<std::vector<double> >& a, std::vector<double>& b, std::vector<double>& x,
		  std::vector<std::vector<double> >& l, std::vector<std::vector<double> >& r, std::vector<double>& y)
{
	std::size_t n = a.size();
	//循环主体(LU)分解
	for(std::size_t i=1; i<=n; ++i)
	{
		//计算R
		for(std::size_t j=i; j<=n; ++j)
		{
			double m1=0;
			for(std::size_t k=1; k<=i-1; ++k)
				m1+=l[i-1][k-1]*r[k-1][j-1];
			r[i-1][j-1]=a[i-1][j-1]-m1;
		}
		//计算Y
		double m2=0;
		for(std::size_t k=1; k<=i-1; ++k)
			m2+=l[i-1][k-1]*y[k-1];
		y[i-1]=b[i-1]-m2;
		//计算L
		for(std::size_t j=i+1; j<=n; ++j)
		{
			double m3=0;
			for(std::size_t k=1; k<=i-1; ++k)
				m3+=l[j-1][k-1]*r[k-1][i-1];
			l[j-1][i-1]=(a[j-1][i-1]-m3)/r[i-1][i-1];
		}
		l[i-1][i-1]=1;  //L对角线元素都变为1
	}
	//计算X
	x[n-1]=y[n-1]/r[n-1][n-1];
	for(int i=n-1; i!=0; i--)
	{	
		double m=0;
		for(std::size_t j=i+1;j<=n;j++)
			m+=r[i-1][j-1]*x[j-1];
		x[i-1]=(y[i-1]-m)/r[i-1][i-1];
	}
}
//-----------------------------------
//判断矩阵对称性 返回值为1对称
bool SymMatrix(std::vector<std::vector<double> >& a)
{
	std::size_t m = a.size();
	std::size_t n = a[0].size();
	if(m!=n) return 0;
	for(std::size_t i=0; i!=n; ++i)
		for(std::size_t j=i; j!=n; ++j)
			if(a[i][j]!=a[j][i]) return 0;
	return 1;
}
//-----------------------------------
//改进平方根法解系数矩阵对称正定的方程组(a增广矩阵）
void SqrtMethod(std::vector<std::vector<double> >& a, std::vector<double>& x)
{
	std::size_t n = a.size();
	std::vector<double> temp1(n,0);
	std::vector<double> temp2(n+1,0);
	std::vector<std::vector<double> > l;
	std::vector<std::vector<double> > r;
	for(std::size_t i=0; i!=n; ++i) l.push_back(temp1);
	for(std::size_t i=0; i!=n; ++i) r.push_back(temp2);
	//循环主体
	for(std::size_t i=1; i<=n; ++i)
	{
		for(std::size_t j=i; j<=n+1; ++j)
		{
			double sum = 0;
			for(std::size_t k=1; k<=i-1; ++k)
				sum+=l[i-1][k-1]*r[k-1][j-1];
			r[i-1][j-1]=a[i-1][j-1]-sum;
		}
		for(std::size_t j=i+1; j<=n; ++j)
			l[j-1][i-1]=r[i-1][j-1]/r[i-1][i-1];
	}
	//计算X
	for(int i=n; i>=1; --i)
	{
		double m=0;
		for(std::size_t k=i+1; k<=n; ++k)
			m+=r[i-1][k-1]*x[k-1];
		x[i-1]=(r[i-1][n]-m)/r[i-1][i-1];
	}
}
//--------------------------------------
//拉格朗日插值法
// x插值节点，ｙ插值节点处的值，ｐ所求点，函数返回ｐ点的值
double LGLR(const std::vector<double>& x, const std::vector<double>& y, double p)
{
	std::size_t n=x.size();
	if(n==0||n!=y.size()) { std::cerr<<"input error!"; return 0; }
	double value=0;
	for(std::size_t i=0; i!=n; ++i)
	{
		double m=y[i];
		for(std::size_t j=0; j!=n; ++j)
			if(i!=j) m*=(p-x[j])/(x[i]-x[j]);
		value+=m;
	}
	return value;
}
//---------------------------------------
//牛顿插值多项式
//N(x)=f[x0]+(x-x0){f[x0,x1]+(x-x1){f[x0,x1,x2]+(x-x2){...{f[x0,...xn-1]+(x-xn-1)f[x0,...xn]}...} 
class CNewton 
{ 
	double *f[2]; 
	double *x; 
	int max; 
	int n; 
public: 
	CNewton(int MaxN);//MaxN 为最大插值点数 可任意设定 
	~CNewton(); 
	void InsertPoint(const Point&); 
	double GetValue(double); 
}; 
CNewton::CNewton(int MaxN) 
{ 
	max=MaxN+1; 
	n=0; 
	x=new double[max]; 
	f[0]=new double[max]; 
	f[1]=new double[max]; 
	assert(x!=NULL); 
	assert(f[0]!=NULL); 
	assert(f[1]!=NULL); 
} 
CNewton::~CNewton() 
{ 
	if(x) delete[]x; 
	if(f[0]) delete[]f[0]; 
	if(f[1]) delete[]f[1]; 
} 
void CNewton::InsertPoint(const Point& p) 
{ 
	double fw; 
	assert(n<max); 
	//重复点检查 
	for(int i=0;i<n;++i) 
	  if(fabs(p.getX()-x[i])<1e-5) std::cout<<"wrong input"; 
	x[n]=p.getX(); 
	fw=p.getY(); 
	for(int i=1;i<=n;++i) 
	{ 
		double tmp=fw; 
		fw=(fw-f[1][i-1])/(x[n]-x[n-i]); 
		f[1][i-1]=tmp; 
	} 
	f[0][n]=f[1][n]=fw; 
	n++; 
} 
double CNewton::GetValue(double X) 
{ 
	if(n==0) return 0.0; 
	double s=f[0][n-1]; 
	for(int i=n-2;i>=0;--i) 
		s=s*(X-x[i])+f[0][i]; 
	return s; 
} 
//------------------------------------
//三次样条函数
// boundary边界条件：１为第一类边界条件，即知道边界上的一阶导；２为第二类边界条件，即知道边界上的二阶导
// 　　　　　 3为周期为xn-x0的周期函数
// s0,sn为边界上的条件，当boundary为２且s0,sn都为０时，又称为自然边界条件
// M二阶导数，SLOP一阶导数
void Spline(const std::vector<Point>& a, std::vector<double>& M, std::vector<double>& slope,
			int boundary, double s0=0, double sn=0)
{
	assert(boundary>0&&boundary<4);
	std::size_t n = a.size();
	assert(n==M.size());

	std::vector<double> h(n-1,0);
	for(std::size_t i=0; i!=n-1; ++i)
		h[i]=a[i+1].getX()-a[i].getX();

	std::vector<double> u(n-2,0);
	std::vector<double> t(u);
	std::vector<double> diagonal(n,2);
	for(std::size_t i=0; i!=n-2; ++i)
	{
		u[i]=h[i]/(h[i]+h[i+1]);
		t[i]=1-u[i];
	}
	
	std::vector<double> d(n-2,0);
	for(std::size_t i=0; i!=n-2; ++i)
		d[i]=6*((a[i+2].getY()-a[i+1].getY())/h[i+1]-(a[i+1].getY()-a[i].getY())/h[i])/(h[i]+h[i+1]);
	
	switch(boundary)
	{
	case 1: //第一类边界条件
		{	
			double D0=6*((a[1].getY()-a[0].getY())/h[0]-s0)/h[0];
			double Dn=6*(sn-(a[n-1].getY()-a[n-2].getY())/h[n-2])/h[n-2];
			d.insert(d.begin(),D0);
			d.push_back(Dn);
			u.push_back(1);
			t.insert(t.begin(),1);
			Chase(u,diagonal,t,d,M);
		}
			for(std::size_t i=0; i!=M.size(); ++i)
				std::cout<<'M'<<i<<": "<<M[i]<<std::endl;
			break;
	case 2: //第二类边界条件
			M.resize(n-2);
			u.erase(u.begin());
			t.pop_back();
			diagonal.resize(n-2);
			Chase(u,diagonal,t,d,M);
			M.insert(M.begin(),s0);
			M.push_back(sn);
			for(std::size_t i=0; i!=M.size(); ++i)
				std::cout<<'M'<<i<<": "<<M[i]<<std::endl;
			break;
	case 3: //第三类边界条件
			M.resize(n-1);
			std::vector<std::vector<double> > matrix;
			std::vector<double> vec(n-1,0);
			for(std::size_t i=0; i!=n-1; ++i)
				matrix.push_back(vec);
			matrix[0][n-1]=u[0];
			matrix[n-1][0]=t[n-3];
			for(std::size_t i=0; i!=n-1; ++i)
			{
				matrix[i][i]=2;
				matrix[i][i+1]=t[i];
				matrix[i+1][i]=u[i];
			}
			double dn=6*((a[0].getY()-a[n-1].getY())/h[0]-(a[n-1].getY()-a[n-2].getY())/h[n-2])/(h[n-2]+h[0]);
			d.push_back(dn);
			Gauss(matrix,d,M);
			for(std::size_t i=0; i!=M.size(); ++i)
				std::cout<<'M'<<i<<": "<<M[i]<<std::endl;
	};
	//计算一阶导
	slope.push_back(-h[0]*(2*M[0]+M[1])/6+(a[1].getY()-a[0].getY())/h[0]);
	for(std::size_t i=1; i!=n; ++i)
		slope.push_back(h[i-1]*(2*M[i]+M[i-1])/6+(a[i].getY()-a[i-1].getY())/h[i-1]);
	for(std::size_t i=0; i!=slope.size(); ++i)
				std::cout<<"slope"<<i<<": "<<slope[i]<<std::endl;
}
//-------------------------------------------
//最小二乘函数逼近
//适用于基函数为x^k的，其中k>=0
//a为离散点，b为每点的权系数，N为拟合的多项式次数，c为所求拟合多项式的系数
void LSM(const std::vector<Point>& a, const std::vector<double>& b, std::size_t N, std::vector<double>& c)
{