椭圆型微分方程

#include<stdio.h>
#include<math.h>
#define m 4
#define n 9
#define e 1.0e-10
#define w 1.25
int Jacobi(int a[n][n],double b[n],double x[n]);//Jacobi迭代法
int SOR(int a[n][n],double b[n],double x[n]); //SOR迭代法函数
int Gauss_Seidel(int a[n][n],double b[n],double x[n]); //Gauss-Seidel迭代法函数
double f(double a,double b);//边值函数
void main()
{
printf("***椭圆型微分方程数值解问题(Dirichlet边值条件)***\n\n");
double x[m+1],y[m+1];
double U[m+1][m+1];
double K[n];
double u;//精确解
double xx1[n],xx2[n],xx3[n];
int a[n][n];//系数矩阵
double h;
h=1.0/m;
int i,j;


for(i=0;i<=m;i++)
{
x[i]=i*h;
y[i]=i*h;
}
for(i=0;i<=m;i++)//Dirichlet边界值
{
U[i][0]=f(x[i],y[0]);
U[i][m]=f(x[i],y[m]);
U[0][i]=f(x[0],y[i]);
U[m][i]=f(x[m],y[i]);
}

for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
a[i][j]=0;
K[i]=0;
}


for(i=0;i<n;i++)//系数矩阵
for(j=0;j<n;j++)
{
if(i==j) a[i][j]=4;
else if(j-i==1&&j%(m-1)!=0)
{
a[i][j]=-1;
a[j][i]=-1;
}
else if(j-i==m-1)
{
a[i][j]=-1;
a[j][i]=-1;
}
}

for(i=0;i<n;i++)//将K矩阵赋值
{

if(i==0) K[0]=U[1][0]+U[0][1];
else if(i>0&&i<m-2) K[i]=U[i+1][0];
else if(i==m-2) K[i]=U[i+1][0]+U[i+2][1];

else if(i==n-m+1) K[i]=U[0][m-1]+U[1][m];
else if(i>n-m+1&&i<n-1) K[i]=U[i-(m-2)*(m-1)+1][m];
else if(i==n-1) K[i]=U[m-1][m]+U[m][m-1];
}
for(j=1;j<m-2;j++)
{
K[(m-1)*j]=U[0][j+1];
K[(m-1)*(j+1)-1]=U[m][j+1];
}


/*
//测试结果正确性与否的已知数据如下(m=4,n=9)：
K[0]=100;
K[1]=20;
K[2]=20;
K[3]=80;
K[4]=0;
K[5]=0;
K[6]=260;
K[7]=180;
K[8]=180;
*/


printf("向量K序列值如下：\n");
for(i=0;i<n;i++)
printf("K[%d]=%-10.2lf\n",i+1,K[i]);


printf("系数矩阵A如下：\n");
for(i=0;i<n;i++)
{

for(j=0;j<n;j++)
{
printf("%-4d",a[i][j]);
}
printf("\n");
}

Jacobi(a,K,xx1);//进行Jacobi迭代
Gauss_Seidel(a,K,xx2);//进行Gauss_Seidel迭代
SOR(a,K,xx3);//进行SOR迭代

printf("迭代求的解为:\n");
printf("U[i] 精确解 Jacobi Gauss_Seidel SOR\n");
for(i=1;i<=n;i++) //打印方程组的解
{
u=f(x[(i-1)%(m-1)+1],y[(i+m-2)/(m-1)]);//精确解
printf("U[%-2d] %-15.10lf %-15.10lf %-15.10lf %-15.10lf\n",i,u,xx1[i-1],xx2[i-1],xx3[i-1]);
}


/*
printf("\n");
printf("迭代求法误差为:\n");
printf("U[i] 精确解 Jacobi Gauss_Seidel SOR\n");
for(i=1;i<=n;i++) //打印误差
{
u=f(x[(i-1)%(m-1)+1],y[(i+m-2)/(m-1)]);
printf("u[%-2d] %e %e %e %-15.10lf\n",i,fabs(u-u),fabs(u-xx1[i-1]),fabs(u-xx2[i-1]),fabs(u-xx3[i-1]));
}
printf("\n");
*/


}
int Jacobi(int a[n][n]

,double b[n],double x[n])
{
int i,j;
double Fan8;
//定义向量无穷范数
double sum1,sum2;//定义累加器
int count=0;//计数器
double x0[n],x1[n];
for(i=0;i<n;i++)
x0[i]=1;//读入初值序列
for(i=0;i<n;i++)
x1[i]=x0[i];//初始化

do
{
for(i=0;i<n;i++)
x0[i]=x1[i];
for(i=1;i<=n;i++)
{
for(j=1,sum1=0;j<=i-1;j++)
sum1=sum1+a[i-1][j-1]*x0[j-1];
for(j=i+1,sum2=0;j<=n;j++)
sum2=sum2+a[i-1][j-1]*x0[j-1];
x1[i-1]=(b[i-1]-sum1-sum2)/a[i-1][i-1];//Jacobi迭代公式
}
count++;
for(Fan8=fabs(x1[0]-x0[0]),i=0;i<n;i++)
if(Fan8<fabs(x1[i]-x0[i]))
Fan8=fabs(x1[i]-x0[i]);//进行无穷范数运算

}while(Fan8>e);

for(i=0;i<n;i++)
x[i]=x1[i];//将结果返回
printf("Jacobi方法迭代次数为 %d次!\n\n",count);
return 0;


}


int SOR(int a[n][n],double b[n],double x[n])
{
int i,j;
int count=0;//计数器
double Fan8;//定义向量无穷范数
double sum1,sum2;//定义累加器
double x0[n],x1[n];
//初值序列
for(i=0;i<n;i++)
x0[i]=1;
for(i=0;i<n;i++)
x1[i]=x0[i];//初始化
do
{
for(i=0;i<n;i++)
x0[i]=x1[i];
for(i=1;i<=n;i++)
{
for(j=1,sum1=0;j<=i-1;j++)
sum1=sum1+a[i-1][j-1]*x1[j-1];
for(j=i,sum2=0;j<=n;j++)
sum2=sum2+a[i-1][j-1]*x0[j-1];
x1[i-1]=x0[i-1]+w*(b[i-1]-sum1-sum2)/a[i-1][i-1];//SOR迭代公式
}
count++;
for(Fan8=fabs(x1[0]-x0[0]),i=0;i<n;i++)
if(Fan8<fabs(x1[i]-x0[i]))
Fan8=fabs(x1[i]-x0[i]);//进行无穷范数运算
}while(Fan8>e);
for(i=0;i<n;i++)
x[i]=x1[i];//将结果返回
printf("SOR方法迭代次数为 %d次!\n\n",count);
return 0;
}

int Gauss_Seidel(int a[n][n],double b[n],double x[n])
{
int i,j;
int count=0;
double Fan8;//定义向量无穷范数
double sum1,sum2;//定义累加器
double x0[n],x1[n];
//初值序列
for(i=0;i<n;i++)
x0[i]=1;
for(i=0;i<n;i++)
x1[i]=x0[i];//初始化

do
{
for(i=0;i<n;i++)
x0[i]=x1[i];
for(i=1;i<=n;i++)
{
for(j=1,sum1=0;j<=i-1;j++)
sum1=sum1+a[i-1][j-1]*x1[j-1];
for(j=i+1,sum2=0;j<=n;j++)
sum2=sum2+a[i-1][j-1]*x0[j-1];
x1[i-1]=(b[i-1]-sum1-sum2)/a[i-1][i-1];//Gauss-Seidel迭代公式
}
count++;
for(Fan8=fabs(x1[0]-x0[0]),i=0;i<n;i++)
if(Fan8<fabs(x1[i]-x0[i]))
Fan8=fabs(x1[i]-x0[i]);//进行无穷范数运算
}while(Fan8>e);
for(i=0;i<n;i++)
x[i]=x1[i];//将结果返回
printf("Gauss_Seidel方法迭代次数为%d次!\n\n",count);
return 0;
}


double f(double a,double b)//边值函数
{
double fun;
fun=pow(a,2)-pow(b,2);
return fun;
}