📄 jifen.c
字号:
/*jifen.c 包括复化梯形积分法、复化Simpson积分法、Gauss_Legendre积分法和Guass_Chebyshev积分法*/
#include "stdio.h"
#include "math.h"
#define MAX 5000 /**控制矩阵最大维数**/
#define eps 1.0e-10 /****控制误差要求****/
#define pi 3.141592653589793
double fftx();
double simp();
double Gauss_Legendre();
double Guass_Chebyshev();
double g(double);
double Liezhu_Guass(int);
double A[MAX][MAX],B[MAX],u[MAX];
main()
{
fftx();
simp();
Gauss_Legendre();
Guass_Chebyshev();
}
/*****************计算函数g(x)****************/
double g(double x)
{
return exp(4*x)+(exp(x+4)-exp(-x-4))/(x+4);
}
/********************************************/
/*****************复化梯形积分法*************/
double fftx()
{
FILE *fp;
int i,j,n;
double x,y,h,ux;
double e=0.0;
n=1;
while(n<=MAX) /*计算方程组的系数矩阵和右端向量B*/
{
h=2.0/n;
for(i=0;i<=n;i++)
{
y=-1+i*h;
for(j=1;j<n;j++)
{
x=-1+j*h;
A[i][j]=2*exp(x*y);
}
A[i][0]=exp(-y);
A[i][n]=exp(y);
}
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
A[i][j]=A[i][j]*h/2.0;
for(i=0;i<=n;i++)
{
B[i]=g(-1+i*h);
A[i][i]+=1;
}
Liezhu_Guass(n); /*解方程组求解数值解*/
e+=pow((exp(4*(-1))-u[0]),2);/*计算误差*/
e+=pow((exp(4*1)-u[n]),2);
for(i=1;i<n;i++)
{
ux=exp(4*(-1+i*h));
e+=2.0*pow((ux-u[i]),2);
}
e=e*h/2.0;
if(e<=eps)
break;
else
n=n*2;
}
if((fp=fopen("C:\\\\Documents and Settings\\\\MN\\\\桌面\\\\jifen\\\\x1.txt","w"))==NULL)
{
printf("Can not open out.txt\n");
exit(1);
}
for(i=0;i<=n;i++)
{
x=-1+h*i;
fprintf(fp,"%lf ",x);
}
fclose(fp);
if((fp=fopen("C:\\\\Documents and Settings\\\\MN\\\\桌面\\\\jifen\\\\u1.txt","w"))==NULL)
{
printf("Can not open out.txt\n");
exit(1);
}
for(i=0;i<=n; i++)
fprintf(fp,"%lf ",u[i]);
fclose(fp);
return u[n];
}
/*******************************************************************/
/****************************复化Simpson积分法**********************/
double simp()
{
FILE *fp;
int i,j,m;
double x,y,h,ux;
double e=0.0;
m=1;
while(2*m+1<=MAX) /*计算方程组的系数矩阵和右端向量B*/
{
h=2.0/(2*m);
for(i=0;i<=2*m;i++)
{
y=-1+i*h;
for(j=1;j<m;j++)
{
x=-1+(2*j-1)*h;
A[i][2*j-1]=exp(x*y)*4;
x=-1+(2*j)*h;
A[i][2*j]=exp(x*y)*2;
}
A[i][0]=exp(-y);
A[i][2*m]=exp(y);
A[i][2*m-1]=exp((-1+(2*m-1)*h)*y)*4;
}
for(i=0;i<=2*m;i++)
for(j=0;j<=2*m;j++)
A[i][j]=A[i][j]*h/3.0;
for(i=0;i<=2*m;i++)
{
B[i]=g(-1+i*h);
A[i][i]+=1;
}
Liezhu_Guass(2*m); /*解方程组求解数值解*/
e+=pow((exp(4*(-1))-u[0]),2); /*计算误差*/
e+=pow((exp(4*1)-u[2*m]),2);
for(i=1;i<=m;i++)
{
ux=exp(4*(-1+(2*i-1)*h));
e+=4.0*pow((ux-u[2*i-1]),2);
}
for(i=1;i<m;i++)
{
ux=exp(4*(-1+(2*i)*h));
e+=2.0*pow((ux-u[2*i]),2);
}
e=e*h/3.0;
if(e<=eps)
break;
else
m=m*2;
}
if((fp=fopen("C:\\\\Documents and Settings\\\\MN\\\\桌面\\\\jifen\\\\x2.txt","w"))==NULL)
{
printf("Can not open out.txt\n");
exit(1);
}
for(i=0;i<=2*m;i++)
{
x=-1+h*i;
fprintf(fp,"%lf ",x);
}
fclose(fp);
if((fp=fopen("C:\\\\Documents and Settings\\\\MN\\\\桌面\\\\jifen\\\\u2.txt","w"))==NULL)
{
printf("Can not open out.txt\n");
exit(1);
}
for(i=0;i<=2*m; i++)
fprintf(fp,"%lf ",u[i]);
fclose(fp);
return u[2*m];
}
/*******************************************************************/
/********************Gauss_Legendre积分法***************************/
double Gauss_Legendre()
{
int i,j,n;
double e=0.0;
double x[9],Ai[9];
FILE *fp;
n=9;
x[0]=-0.9681602395; /*积分节点*/
x[8]=0.9681602395;
x[1]=-0.8360311073;
x[7]=0.8360311073;
x[2]=-0.6133714327;
x[6]=0.6133714327;
x[3]=-0.3242534234;
x[5]=0.3242534234;
x[4]=0.0;
Ai[0]=Ai[8]=0.0812743884; /*积分系数*/
Ai[1]=Ai[7]=0.1806481607;
Ai[2]=Ai[6]=0.2606106964;
Ai[3]=Ai[5]=0.3123470770;
Ai[4]=0.3302393550;
for(i=0;i<n;i++) /*计算系数矩阵A和右端向量B*/
{
for(j=0;j<n;j++)
{
A[i][j]=Ai[j]*exp(x[i]*x[j]);
}
}
for(i=0;i<n;i++)
{
B[i]=g(x[i]);
A[i][i]+=1;
}
Liezhu_Guass(n-1); /*解方程组求解数值解*/
for(i=0;i<n;i++)
e+=Ai[i]*pow((exp(4*x[i])-u[i]),2);/*计算误差*/
if((fp=fopen("C:\\\\Documents and Settings\\\\MN\\\\桌面\\\\jifen\\\\x3.txt","w"))==NULL)
{
printf("Can not open out.txt\n");
exit(1);
}
for(i=0;i<n;i++)
fprintf(fp,"%lf ",x[i]);
fclose(fp);
if((fp=fopen("C:\\\\Documents and Settings\\\\MN\\\\桌面\\\\jifen\\\\u3.txt","w"))==NULL)
{
printf("Can not open out.txt\n");
exit(1);
}
for(i=0;i<n; i++)
fprintf(fp,"%lf ",u[i]);
fclose(fp);
return u[n];
}
/****************************************************************************/
/***********************Guass-Chebyshev积分法********************************/
double Guass_Chebyshev()
{
int i,j,n;
double e=0.0;
double x[MAX];
FILE *fp;
n=1;
while(n<=MAX)
{
for(i=0;i<n;i++) /*计算节点*/
{
x[i]=cos((2*(n-i-1)+1)*pi/(2*n));
}
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
A[i][j]=pi*sqrt(1-pow(x[j],2))*exp(x[i]*x[j])/n;
}
}
for(i=0;i<n;i++)
{
B[i]=g(x[i]);
A[i][i]+=1; /*计算系数矩阵A和右端向量B*/
}
Liezhu_Guass(n-1); /*解方程组求解数值解*/
for(i=0;i<n;i++) /*计算误差*/
{
e+=pow((exp(4*x[i])-u[i]),2);
}
e=e*pi/n;
if(e<=eps)
break;
else
n=n*2;
}
if((fp=fopen("C:\\\\Documents and Settings\\\\MN\\\\桌面\\\\jifen\\\\x4.txt","w"))==NULL)
{
printf("Can not open out.txt\n");
exit(1);
}
for(i=0;i<n;i++)
fprintf(fp,"%lf ",x[i]);
fclose(fp);
if((fp=fopen("C:\\\\Documents and Settings\\\\MN\\\\桌面\\\\jifen\\\\u4.txt","w"))==NULL)
{
printf("Can not open out.txt\n");
exit(1);
}
for(i=0;i<n; i++)
fprintf(fp,"%lf ",u[i]);
fclose(fp);
return u[n];
}
/****************************************************************************/
/**********************列主高斯消元法解方程组********************************/
double Liezhu_Guass(int n)
{
int i,j,k,max;
double t,m;
for(k=0;k<n;k++)
{
max=k;
for(i=k+1;i<=n;i++) /*选绝对值最大的元素的行号*/
if(fabs(A[max][k])<fabs(A[i][k]))
max=i;
for(j=k;j<=n;j++)
{
t=A[k][j]; /*交换A[k][]行和A[max][]行*/
A[k][j]=A[max][j];
A[max][j]=t;
}
t=B[k]; /*交换B[k]行和B[max]行*/
B[k]=B[max];
B[max]=t;
for(i=k+1;i<=n;i++) /*高斯消元，将A化为上三角阵*/
{
m=A[i][k]/A[k][k];
for(j=k;j<=n;j++)
{
A[i][j]-=m*A[k][j];
}
B[i]-=m*B[k];
}
}
u[n]=B[n]/A[n][n]; /*回代*/
for(k=n-1;k>=0;k--)
{
u[k]=0.0;
for(j=k+1;j<=n;j++)
{
u[k]+=A[k][j]*u[j];
}
u[k]=(B[k]-u[k])/A[k][k];
}
return u[n+1];
}
/*********************************************************************/⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -