2次元ラプラス方程式の差分法のプログラム

先週、渡したプリントに書かれていた練習問題の解答を示します。
ソースプログラム

#include 
#define MAXDIM 510

void initialize(double u[][MAXDIM], int flag[][MAXDIM]);
void set_cordinate(int n,double x[][MAXDIM],double y[][MAXDIM]);
void set_boundary_wall_pot(int n, double u[][MAXDIM], int f[][MAXDIM]);
void set_circle(double x0, double y0, double r, double v,
		int n, double x[][MAXDIM], double y[][MAXDIM],
	        double u[][MAXDIM], int f[][MAXDIM]);
void plot_3d(FILE *gp, char *data_file,
	     double rot_x, double rot_z);

/*=====================================================================*/
/*   main                                                              */
/*=====================================================================*/
int main(){
   double u[MAXDIM][MAXDIM];
   double x[MAXDIM][MAXDIM], y[MAXDIM][MAXDIM];
   int flag[MAXDIM][MAXDIM];
   int i,j;
   int k, iteration;
   int nlat;
   char *plot_data_file;
   FILE *fp;
   FILE *gp;
   
   plot_data_file="result.txt";
   
   nlat=100;
   
   initialize(u, flag);
   set_cordinate(nlat,x,y);
   set_boundary_wall_pot(nlat, u, flag);
   set_circle(0.3, 0.3, 0.2,  20.0, nlat, x, y, u, flag);
   set_circle(0.7, 0.6, 0.03, -20.0, nlat, x, y, u, flag);
   
   iteration=2000;
   
   for(k=1; k<=iteration;k++){
     for(j=1; j<=nlat-1; j++){
       for(i=1; i<=nlat-1;i++){
	 if(flag[i][j] != 1){
	   u[i][j]=0.25*(u[i][j-1]+u[i][j+1]+u[i-1][j]+u[i+1][j]);
	 }
       }
     }
   }
   
   
   fp=fopen(plot_data_file,"w");
   
   for(j=0; j<=nlat; j++){
     for(i=0; i<=nlat; i++){
       fprintf(fp, "%f\t%f\t%f\n",x[i][j],y[i][j],u[i][j]);
     }
     fprintf(fp,"\n");
   }
   
   fclose(fp);
   
   
   gp = popen("gnuplot -persist","w");
   plot_3d(gp, plot_data_file, 30, 30);
   plot_3d(gp, plot_data_file, 45, 30);
   plot_3d(gp, plot_data_file, 60, 30);
   plot_3d(gp, plot_data_file, 75, 30);
   pclose(gp);

   return 0;
   
}

/*=====================================================================*/
/*   initialize function                                               */
/*=====================================================================*/
void initialize(double u[][MAXDIM], int flag[][MAXDIM]){

  int i, j;

  for(i=0; i<=MAXDIM-1; i++){
    for(j=0; j<=MAXDIM-1; j++){
      u[i][j]=0.0;
      flag[i][j]=0;
    }
  }
}

/*=====================================================================*/
/*   set coordinate                                                    */
/*=====================================================================*/
void set_cordinate(int n,double x[][MAXDIM],double y[][MAXDIM]){

  int i, j;

  for(i=0; i<=n; i++){
    for(j=0; j<=n; j++){
      x[i][j]=(double)i/n;
      y[i][j]=(double)j/n;
    }
  }
}

/*=====================================================================*/
/*   set external boudary points                                       */
/*=====================================================================*/
void set_boundary_wall_pot(int n, double u[][MAXDIM], int f[][MAXDIM]){

  int i;

  for(i=0; i<=n; i++){
    u[0][i]=0.0;
    f[0][i]=1;
    
    u[n][i]=0.0;
    f[n][i]=1;

    u[i][0]=0.0;
    f[i][0]=1;

    u[i][n]=0.0;
    f[i][n]=1;
  }
  
}

/*=====================================================================*/
/*   set circler boundary                                              */
/*=====================================================================*/
void set_circle(double x0, double y0, double r, double v,
		int n, double x[][MAXDIM], double y[][MAXDIM],
		double u[][MAXDIM], int f[][MAXDIM]){
  
  int i, j;
  
  for(i=0; i<=n; i++){
    for(j=0; j<=n; j++){
      if((x[i][j]-x0)*(x[i][j]-x0)+(y[i][j]-y0)*(y[i][j]-y0) < r*r){
	u[i][j]=v;
	f[i][j]=1;
      }
    }
  }
}


/*==========================================================*/
/*   3D plot                                                */
/*==========================================================*/
void plot_3d(FILE *gp, char *data_file,
	     double rot_x, double rot_z){
  
  
  /* == flowing lines make a graph by using gnuplot == */
  
  printf("start plot 3D\n");
  fprintf(gp, "reset\n");
  fprintf(gp, "set terminal postscript eps color\n");
  fprintf(gp, "set output \"graph.eps\"\n");
  fprintf(gp, "set view %f,%f\n", rot_x, rot_z);
  fprintf(gp, "set contour base\n");
  fprintf(gp, "set cntrparam levels 20\n");
  fprintf(gp, "set hidden3d\n");
  fprintf(gp, "set nokey\n");
  fprintf(gp, "splot \"%s\" with lines\n", data_file);
  fprintf(gp, "set terminal x11\n");
  fprintf(gp, "replot\n");
  
}
結果の図

**図 1:**最初の図
**図 2:**2番目の図
**図 3:**3番目の図
**図 4:**4番目の図
last update:
2次元ラプラス方程式の差分法のプログラム

目次

ソースプログラム

結果の図
