nownow.cpp

分别采用复化梯形公式

#include "stdio.h"
#include "math.h"

/********************************复化梯形积分法*****************************************/
void jifentixing()
{
double a=-1.0;          /*定义a,b积分区间*/
double b=1.0;
double middle,hm,error; /*定义h,误差error*/
int n=2000;             /*定义节点数*/
int i,j,k,im,tt;
#define len 2000      
double x[len+1],g[len+1],u[len+1];/*定义x,g(x)和u(x)*/
double aa[len+1][len+1];          /*用于求解方程组*/
double a3[len+2][len+2],l1[len+2][len+2],u1[len+2][len+2];/*用于求解方程组*/
double u0[len+2],y[len+2],l[len+2];/*用于求解方程组*/
hm=(b-a)/len;/*求解h*/
for (i=0;i<n+1;i++)
{ x[i]=a+i*hm;}/*求解xi*/

for (i=0;i<n+1;i++)
{ g[i]=exp(4*x[i])+(exp(x[i]+4)-exp(-(x[i]+4)))/(x[i]+4);}/*求解g(xi)*/



for (i=0;i<n+1;i++)
	{
		aa[i][0]=hm/2.0*exp(x[i]*a);
        if (i==0)
		{aa[i][0]=aa[i][0]+1;}
	} 
for (j=1;j<n;j++)
{ 
	for (i=0;i<n+1;i++)
	{
		aa[i][j]=hm*exp(x[i]*x[j]);     /*求解方程组系数矩阵*/
        if (i==j)
		{aa[i][j]=aa[i][j]+1;}
	}
}
for (i=0;i<n+1;i++)
	{
		aa[i][n]=hm/2.0*exp(x[i]*b);
        if (i==n)
		{aa[i][j]=aa[i][j]+1;}
	}
for (i=1;i<n+2;i++)
{
	for (j=1;j<n+2;j++)
	{ a3[i][j]=aa[i-1][j-1];}
}



/***********************利用Doolittle分解法求解方程组*************************************/
for (j=1;j<n+2;j++)
{u1[1][j]=a3[1][j];}

for (i=2;i<n+2;i++)
{l1[i][1]=a3[i][1]/u1[1][1];}

for (im=1;im<n+2;im++)
{
	for (j=im;j<n+2;j++)
	{	middle=0.0;
		for (tt=1;tt<im;tt++)
		{ middle=middle+l1[im][tt]*u1[tt][j];}
		u1[im][j]=a3[im][j]-middle;
	}

	if (im<n+1)
	{
		for (i=im+1;i<n+2;i++)
		{   middle=0.0;
			for (tt=1;tt<im;tt++)
			{middle=middle+l1[i][tt]*u1[tt][im];}
			l1[i][im]=(a3[i][im]-middle)/u1[im][im];
		}
	}
} /*得到分解后的L和U矩阵 */
/*************分解过程结束************************************************/
/*****求解过程*****/
for (i=1;i<n+2;i++)
{ y[i]=g[i-1];}

l[1]=y[1];
for (i=2;i<n+2;i++)
{	middle=0.0;
	for (tt=1;tt<i;tt++)
	{middle=middle+l1[i][tt]*l[tt];}
	l[i]=y[i]-middle;}

u0[n+1]=l[n+1]/u1[n+1][n+1];
for (i=n;i>0;i--)
{	middle=0;
	for (tt=i+1;tt<n+2;tt++)
	{middle=middle+u1[i][tt]*u0[tt];}
	u0[i]=(l[i]-middle)/u1[i][i];
}

for (i=0;i<n+1;i++)
{ u[i]=u0[i+1];}
/*求解方程组得到方程组的解ui*/
/*******************************判断误差是否满足条件*******************/
middle=0.0;
for (k=1;k<n;k++)
{middle=middle+2.0*(exp(4*x[k])-u[k])*(exp(4*x[k])-u[k]);}

error=hm/2.0*((exp(4*a)-u[0])*(exp(4*a)-u[0])+(exp(4*b)-u[n])*(exp(4*b)-u[n])+middle);
printf("%13e  ",error);   /*输出结果*/

}

/**************************复化simpson积分方法*************************/
void jifensimpson()
{

double a=-1.0;   /*定义a,b积分区间*/
double b=1.0;
double middle,hm,error;  /*定义h,误差error*/
int m,n;
int i,j,k,im,tt; 
m=40;
n=80;                     /*定义节点数*/
#define len 80
double x[len+1],g[len+1],u[len+1];  /*定义x,g(x)和u(x)*/
double aa[len+1][len+1];                 /*用于求解方程组*/
double a3[len+2][len+2],l1[len+2][len+2],u1[len+2][len+2];   /*用于求解方程组*/
double u0[len+2],y[len+2],l[len+2];      /*用于求解方程组*/
hm=(b-a)/len;     /*求解h*/
for (i=0;i<n+1;i++)
{ x[i]=a+i*hm;}       /*求解xi*/

for (i=0;i<n+1;i++)
{ g[i]=exp(4*x[i])+(exp(x[i]+4)-exp(-(x[i]+4)))/(x[i]+4);}/*求解g(xi)*/



for (i=0;i<n+1;i++)
	{
		aa[i][0]=hm/3.0*exp(x[i]*a);
        if (i==0)
		{aa[i][0]=aa[i][0]+1;}            /*求解方程组系数矩阵*/
	} 
for (j=1;j<n;j=j+2)
{ 
	for (i=0;i<n+1;i++)
	{
		aa[i][j]=(4.0/3.0)*hm*exp(x[i]*x[j]);
        if (i==j)
		{aa[i][j]=aa[i][j]+1;}
	}
}
for (j=2;j<n-1;j=j+2)
{ 
	for (i=0;i<n+1;i++)
	{
		aa[i][j]=(2.0/3.0)*hm*exp(x[i]*x[j]);
        if (i==j)
		{aa[i][j]=aa[i][j]+1;}
	}
}


for (i=0;i<n+1;i++)
	{
		aa[i][n]=hm/3.0*exp(x[i]*b);
        if (i==n)
		{aa[i][j]=aa[i][j]+1;}
	}


for (i=1;i<n+2;i++)
{
	for (j=1;j<n+2;j++)
	{ a3[i][j]=aa[i-1][j-1];}
	   
}


/***********************利用Doolittle分解法求解方程组*************************************/
{u1[1][j]=a3[1][j];}

for (i=2;i<n+2;i++)
{l1[i][1]=a3[i][1]/u1[1][1];}

for (im=1;im<n+2;im++)
{
	for (j=im;j<n+2;j++)
	{	middle=0.0;
		for (tt=1;tt<im;tt++)
		{ middle=middle+l1[im][tt]*u1[tt][j];}
		u1[im][j]=a3[im][j]-middle;
	}

	if (im<n+1)
	{
		for (i=im+1;i<n+2;i++)
		{   middle=0.0;
			for (tt=1;tt<im;tt++)
			{middle=middle+l1[i][tt]*u1[tt][im];}
			l1[i][im]=(a3[i][im]-middle)/u1[im][im];
		}
	}
} /*得到分解后的L和U矩阵 */
/*************分解过程结束************************************************/
/*****求解过程*****/
for (i=1;i<n+2;i++)
{ y[i]=g[i-1];}

l[1]=y[1];
for (i=2;i<n+2;i++)
{	middle=0.0;
	for (tt=1;tt<i;tt++)
	{middle=middle+l1[i][tt]*l[tt];}
	l[i]=y[i]-middle;}

u0[n+1]=l[n+1]/u1[n+1][n+1];
for (i=n;i>0;i--)
{	middle=0;
	for (tt=i+1;tt<n+2;tt++)
	{middle=middle+u1[i][tt]*u0[tt];}
	u0[i]=(l[i]-middle)/u1[i][i];
}

for (i=0;i<n+1;i++)
{ u[i]=u0[i+1];}
/*求解方程组得到方程组的解ui*/
/*******************************判断误差是否满足条件*******************/
middle=0.0;
for (k=1;k<m+1;k++)
{middle=middle+4*(exp(4*x[2*k-1])-u[2*k-1])*(exp(4*x[2*k-1])-u[2*k-1]);}

for (k=1;k<m;k++)
{middle=middle+2*(exp(4*x[2*k])-u[2*k])*(exp(4*x[2*k])-u[2*k]);}

error=(hm/3)*((exp(4*a)-u[0])*(exp(4*a)-u[0])+(exp(4*b)-u[n])*(exp(4*b)-u[n])+middle);
printf("%13e  ",error);/*输出结果*/


}
/*****************************Gauss积分法**************************/
void jifengauss()
{double a=-1.0;     /*定义a,b积分区间*/
double b=1.0;
double middle,error;   /*定义误差error*/
int n;
int i,j,k,im,tt;
double x[8],s[8],ai[8];    /*定义xi,si和ai*/
#define len 7       /*定义节点数*/
double g[len+1],u[len+1];   /*定义g(x)和u(x)*/
double aa[len+1][len+1];    /*用于求解方程组*/
double a3[len+1][len+1],l1[len+1][len+1],u1[len+1][len+1]; /*用于求解方程组*/
double u0[len+1],y[len+1],l[len+1]; /*用于求解方程组*/

x[1]=0.9491079123;
x[2]=-0.9491079123;
x[3]=0.7415311856;/*查表得xi*/
x[4]=-0.7415311856;
x[5]=0.4058451514;
x[6]=-0.4058451514;
x[7]=0;
n=7;
for (i=1;i<n+1;i++)
{s[i]=2/(1-x[i]*x[i]);}  /*求得si*/

ai[1]=0.1294849662;
ai[2]=0.1294849662;
ai[3]=0.2797053915;/*查表得ai*/
ai[4]=0.2797053915;
ai[5]=0.3818300505;
ai[6]=0.3818300505;
ai[7]=0.4179591837;

for (i=1;i<n+1;i++)
{ g[i]=exp(4*x[i])+(exp(x[i]+4)-exp(-(x[i]+4)))/(x[i]+4);} /*求解g(xi)*/


for (j=1;j<n+1;j++)
{ 
	for (i=1;i<n+1;i++)
	{
		aa[i][j]=ai[j]*exp(x[i]*x[j]);     /*求解方程组系数矩阵*/
        if (i==j)
		{aa[i][j]=aa[i][j]+1;}
	}
}

for (i=1;i<n+1;i++)
{
	for (j=1;j<n+1;j++)
	{ a3[i][j]=aa[i][j];}
}



/***********************利用Doolittle分解法求解方程组*************************************/
for (j=1;j<n+1;j++)
{u1[1][j]=a3[1][j];}

for (i=2;i<n+1;i++)
{l1[i][1]=a3[i][1]/u1[1][1];}

for (im=1;im<n+1;im++)
{
	for (j=im;j<n+1;j++)
	{	middle=0.0;
		for (tt=1;tt<im;tt++)
		{ middle=middle+l1[im][tt]*u1[tt][j];}
		u1[im][j]=a3[im][j]-middle;
	}

	if (im<n)
	{
		for (i=im+1;i<n+1;i++)
		{   middle=0.0;
			for (tt=1;tt<im;tt++)
			{middle=middle+l1[i][tt]*u1[tt][im];}
			l1[i][im]=(a3[i][im]-middle)/u1[im][im];
		}
	}
} /*得到分解后的L和U矩阵 */
/*************分解过程结束************************************************/
/*****求解过程*****/
for (i=1;i<n+1;i++)
{ y[i]=g[i];}

l[1]=y[1];
for (i=2;i<n+1;i++)
{	middle=0.0;
	for (tt=1;tt<i;tt++)
	{middle=middle+l1[i][tt]*l[tt];}
	l[i]=y[i]-middle;}

u0[n]=l[n]/u1[n][n];
for (i=n-1;i>0;i--)
{	middle=0;
	for (tt=i+1;tt<n+1;tt++)
	{middle=middle+u1[i][tt]*u0[tt];}
	u0[i]=(l[i]-middle)/u1[i][i];
}

for (i=1;i<n+1;i++)
{ u[i]=u0[i];}
/*求解方程组得到方程组的解ui*/
/*******************************判断误差是否满足条件*******************/
middle=0.0;
for (k=1;k<n+1;k++)
{middle=middle+ai[k]*(exp(4*x[k])-u[k])*(exp(4*x[k])-u[k]);}
error=middle;
printf("%13e  ",error);/*输出结果*/
}


/*********************************主程序************************************/
void main()
{ jifentixing(); /*运行复化梯形积分法*/
jifensimpson();/*运行复化simpson积分法*/
jifengauss();/*运行Gauss积分法*/
}