seidel迭代.cpp

二维稳态导热问题。设材料物性均为常数

#include<iostream>
#include<math.h>
#include<stdio.h>
#define NM 300                 //对称的下半部分未知数个数
void main()
{
	double a[NM][NM],t[NM],s;  //a[NM][NM]为系数矩阵,t[NM]存放对称的下半部分温度
 	double T[450],e=100.0,e1;  //T[]存放所有结点温度
	float dx;                  //dx为步长  
	int N,i,j,k,Itry=0;        //N为x方向结点个数，Itry为迭代次数
	FILE *outfile;
	printf("输入步长(dx/L):\n");
	scanf("%f",&dx);
	N=int(1.0/dx)+2;           //计算x方向结点个数
	printf("N=%d",N);
	for(i=0;i<N*(N+1)/2;i++)   //系数矩阵赋初值
		for(j=0;j<N*(N+1)/2;j++)
			a[i][j]=0;
    for(i=0;i<N*(N+3)/2;i++)   //假设初始温度
		t[i]=0;
	for(i=N;i<N*(N+3)/2;i=i+N) //左边界温度
	    t[i]=1.0;
	for(j=1;j<(N+1)/2;j++)     //内部结点的系数
		for(i=1;i<N-1;i++)
		{
			k=j*N+i;
			a[k][k-1]=1.0;
			a[k][k+1]=1.0;
			a[k][k-N]=1.0;
			a[k][k+N]=1.0;
			a[k][k]=-4;
		}
	while(e>0.0000005)
		{
			e=0.0;
		    for(i=N+1;i<N*(N+1)/2-1;i++)//迭代内部结点的温度
			{ 
			    if(i%N==0||(i+1)%N==0)  continue; //除去边界结点
		     	s=t[i];                           //把迭代前温度赋给s,以便求误差
		     	t[i]=-1/a[i][i]*(a[i][i-1]*t[i-1]+a[i][i+1]*t[i+1]+a[i][i-N]*t[i-N]+a[i][i+N]*t[i+N]);//内部结点的差分方程
				if(i>N*(N-3)/2&&i<((N-1)*N/2-1))  //把与对称轴相临的那一行温度赋给上边界,作为边界条件
				{
					for(j=N*(N+1)/2+1;j<N*(N+3)/2-1;j++)
					t[j]=t[j-2*N];
				}
		    	e1=fabs(s-t[i]);                  //计算最大误差
		    	if(e1>e)  e=e1;
			}
			Itry++;                               //计算迭代次数
		}
	printf("\nItry=%d\n",Itry);
    for(i=0;i<(N+1)*N/2;i++)                       //根据对称性,把所有结点温度写入T[]
			T[i]=t[i];
	for(i=(N+1)*N/2;i<N*N;i++)
	{
		k=2*(i/N-(N-1)/2);
		j=i-k*N;
		T[i]=t[j];
	}
	outfile=fopen("Seidel迭代.txt","w");           //输出到文本文件
	if(outfile==NULL)
	{
		printf("打不开数据文件:分析解.txt!");
		exit(1);
	}
    fprintf(outfile,"步长=%f\n",dx);
	fprintf(outfile,"迭代次数=%d\n",Itry);
	for(j=0;j<N;j++)
		for(i=0;i<N;i++)
		{
			k=j*N+i;
			printf("%1.4f ",T[k]);
			fprintf(outfile,"%1.4f  ",T[k]);
			if((k+1)%N==0)
			{
				printf("\n");
				fprintf(outfile,"\n");
			}
		}
    fclose(outfile);
    


}