simpson.cpp

本人的计算方法课程作业

//simpson.cpp--用复化SIMPSON公式求积分方程的离散解，
//用三次样条函数近似，并与真值比较
#include<iostream>
#include<cmath>
using namespace std;
const int N=8;
//f()定义积分方程左边的多项式
float f(float x)
{
	float g;
	g=(4*x*x*x+5*x*x-2*x+5)/(8*(x+1)*(x+1));
	return g;
}
//max_row()选列主元
int max_row(float a[][N+1],int k)
{
	int i,r;
	float max=fabs(a[k][k]);
	r=k;
	for(i=k+1;i<N+1;i++)
	{
		if(fabs(a[i][k])>max)
		{max=fabs(a[i][k]);
		r=i;}
	}
	return r;
}
//change()交换列元素位置
void change(float a[N+1][N+1],int r,int k)
{
	int j;
	float t;
	for(j=k;j<N+1;j++)
	{
		t=a[k][j];
		a[k][j]=a[r][j];
		a[r][j]=t;
	}
}
//print_err()用于输出错误信息
void print_err()
{
	cout<<"Failure!"<<endl;
}
//guass()按列主元素法解N+1阶线性方程组的解
//系数矩阵存放在二维数组a[][]中，右端向量存放在数组b[]中
void gauss(float a[N+1][N+1],float b[N+1])
{
	int i,j,k;
	float t;
	for(k=0;k<N;k++)
	{
		int r=max_row(a,k);
		if(fabs(a[r][k])<1e-5)
		{print_err();return;}
	    if(k!=r)
		{change(a,r,k);
		t=b[k];b[k]=b[r];b[r]=t;}
		for(i=k+1;i<N+1;i++)                 //消元过程
		{
			for(j=k+1;j<N+1;j++)
				a[i][j]=a[i][j]-a[i][k]*a[k][j]/a[k][k];
			b[i]=b[i]-a[i][k]*b[k]/a[k][k];
		}
	}
	b[N]=b[N]/a[N][N];                       //消元结束，开始回代             
	for(i=N-1;i>=0;i--)
	{	
		float t=0;
		for(j=i+1;j<N+1;j++)
			t+=b[j]*a[i][j];
	    b[i]=(b[i]-t)/a[i][i];
	}
}
//FINDK()找出给定插值数据点所在的区间(x[k-1],x[k]),结果返回下标值k
int FINDK(float xx,float x[])
{
	int k=1,i;
	for(i=1;i<=N;i++)
	{
		if (xx<=x[i])
		{ k=i; goto end;}
	    else
		{ k=i+1;}
	}
end:    ;
	  return k;
 }
//EVASPLINE()用三次样条插值函数近似y(x)的离散解
float EVASPLINE(float m[],float x[],float y[],float xx)
{
	int k;
	float h,S;
	k=FINDK(xx,x);
	h=x[k]-x[k-1];
    if(k>=1)
	   S=(m[k-1]*pow((x[k]-xx),3)/6+m[k]*pow((xx-x[k-1]),3)/6
			+(y[k-1]-m[k-1]*h*h/6)*(x[k]-xx)+
				(y[k]-m[k]*h*h/6)*(xx-x[k-1]))/h;
	return S;
}

//SPLINEM()调用TSS算法计算三次样条插值的参数Mi,系数存放于数组a,b,c中，
//d中是右端向量，计算结果存放于数组m[]中,其它各参数根据第三边界条件确定
float SPLINEM(float x[],float y[],float t)
{
	int i,k;
	float h[N+1],a[N+1],b[N+1],c[N+1],d[N+1],m[N+1];
	float s;
	for(i=0;i<=N;i++)
		m[i]=y[i];
	for(k=1;k<=2;k++)
	{
		for (i=N;i>=k;i--)
			m[i]=(m[i]-m[i-1])/(x[i]-x[i-k]);
	}
	h[1]=x[1]-x[0];
	for (i=1;i<=N-1;i++)
	{
		h[i+1]=x[i+1]-x[i];
	    c[i]=h[i+1]/(h[i]+h[i+1]);
	    a[i]=1-c[i];
	    b[i]=2;
	    d[i]=6*m[i+1];
	}
    d[0]=(-12)*h[1]*(m[3]-m[2])/(x[3]-x[0]);
	d[N]=12*h[N]*(m[N]-m[N-1])/(x[N]-x[N-3]);
	c[0]=-2;
	b[0]=2;
	a[N]=-2;
	b[N]=2;
	m[0]=d[0];
	for(k=2;k<=N+1;k++)
	{
		a[k-1]=a[k-1]/b[k-2];
		b[k-1]=b[k-1]-a[k-1]*c[k-2];
		m[k-1]=d[k-1]-a[k-1]*m[k-2];
	}
	m[N]/=b[N];
	for(k=N-1;k>=0;k--)
	    m[k]=(m[k]-c[k]*m[k+1])/b[k]; 
	s=EVASPLINE(m,x,y,t);
	return s;
}
//realf()用于定义真解函数
float realf(float x)
{
    float t;
	t=1/((x+1)*(x+1));
	return t;
}
void main()
{ 
   float a[N+1][N+1],x[N+1],y[N+1];    //定义存放系数矩阵和右端向量的数组
   int i,k;
   float h=2.0/N;                      //确定步长
   for(i=0;i<=N;i++)
     {
	   x[i]=i*1.0/N;                   //采集数据节点，间距为1/8
       y[i]=-f(x[i]);
     }
   for(i=0;i<=N;i++)                   //由复化simpson公式得到系数矩阵
   {
	   a[i][0]=h/6*(1/(1+x[0])-x[i]);
       a[i][N]=h/6*(1/(1+x[N])-x[i]);
	   for(k=1;k<=N-1;k++)
		   {
			   if(k%2==0)
			       a[i][k]=2*h/6*(1/(1+x[k])-x[i]);
		       else
			       a[i][k]=4*h/6*(1/(1+x[k])-x[i]);
		   }
	   a[i][i]-=1;
   }
   gauss(a,y);
   cout<<"此积分方程的离散近似解为:"<<endl;
   for(i=0;i<=N;i++)                        
   cout<<"("<<x[i]<<" , "<<y[i]<<")"<<endl;
   cout<<"在区间上选取一些点，观察并比较三次样条函数与真值函数的计算结果:"<<endl;
   for(i=1;i<=20;i=i+2)
   {
	   float xi,yi,m,ei;
	   xi=1.0*i/20;
	   m=SPLINEM(x,y,xi);
	   yi=realf(xi);
	   ei=fabs(yi-m);
	   cout<<"SPLINEM("<<xi<<")="<<m;
	   cout<<"   realf("<<xi<<")="<<yi;
	   cout<<"   err="<<ei<<endl;
   }
}