ラプラス方程式を解くプログラム

---

ラプラス方程式

001: #include <stdio.h>
002: #define MAXDIM 510
003: 
004: void initialize(double u[][MAXDIM], int flag[][MAXDIM]);
005: void set_cordinate(int n,double x[][MAXDIM],double y[][MAXDIM]);
006: void set_boundary_wall_pot(int n, double u[][MAXDIM], int f[][MAXDIM]);
007: void set_circle(double x0, double y0, double r, double v,
008:         int n, double x[][MAXDIM], double y[][MAXDIM],
009:             double u[][MAXDIM], int f[][MAXDIM]);
010: void plot_3d(FILE *gp, char *data_file,
011:          double rot_x, double rot_z);
012: 
013: /*=====================================================================*/
014: /*   main                                                              */
015: /*=====================================================================*/
016: int main(){
017:    double u[MAXDIM][MAXDIM];
018:    double x[MAXDIM][MAXDIM], y[MAXDIM][MAXDIM];
019:    int flag[MAXDIM][MAXDIM];
020:    int i,j;
021:    int k, iteration;
022:    int nlat;
023:    char *plot_data_file;
024:    FILE *fp;
025:    FILE *gp;
026:    
027:    plot_data_file="result.txt";
028:    
029:    nlat=100;
030:    
031:    initialize(u, flag);
032:    set_cordinate(nlat,x,y);
033:    set_boundary_wall_pot(nlat, u, flag);
034:    set_circle(0.3, 0.3, 0.2,  20.0, nlat, x, y, u, flag);
035:    set_circle(0.7, 0.6, 0.03, -20.0, nlat, x, y, u, flag);
036:    
037:    iteration=2000;
038:    
039:    for(k=1; k<=iteration;k++){
040:      for(j=1; j<=nlat-1; j++){
041:        for(i=1; i<=nlat-1;i++){
042:      if(flag[i][j] != 1){
043:        u[i][j]=0.25*(u[i][j-1]+u[i][j+1]+u[i-1][j]+u[i+1][j]);
044:      }
045:        }
046:      }
047:    }
048:    
049:    
050:    fp=fopen(plot_data_file,"w");
051:    
052:    for(j=0; j<=nlat; j++){
053:      for(i=0; i<=nlat; i++){
054:        fprintf(fp, "%f\t%f\t%f\n",x[i][j],y[i][j],u[i][j]);
055:      }
056:      fprintf(fp,"\n");
057:    }
058:    
059:    fclose(fp);
060:    
061:    
062:    gp = popen("gnuplot -persist","w");
063:    plot_3d(gp, plot_data_file, 30, 30);
064:    plot_3d(gp, plot_data_file, 45, 30);
065:    plot_3d(gp, plot_data_file, 60, 30);
066:    plot_3d(gp, plot_data_file, 75, 30);
067:    pclose(gp);
068: 
069:    return 0;
070:    
071: }
072: 
073: /*=====================================================================*/
074: /*   initialize function                                               */
075: /*=====================================================================*/
076: void initialize(double u[][MAXDIM], int flag[][MAXDIM]){
077: 
078:   int i, j;
079: 
080:   for(i=0; i<=MAXDIM-1; i++){
081:     for(j=0; j<=MAXDIM-1; j++){
082:       u[i][j]=0.0;
083:       flag[i][j]=0;
084:     }
085:   }
086: }
087: 
088: /*=====================================================================*/
089: /*   set coordinate                                                    */
090: /*=====================================================================*/
091: void set_cordinate(int n,double x[][MAXDIM],double y[][MAXDIM]){
092: 
093:   int i, j;
094: 
095:   for(i=0; i<=n; i++){
096:     for(j=0; j<=n; j++){
097:       x[i][j]=(double)i/n;
098:       y[i][j]=(double)j/n;
099:     }
100:   }
101: }
102: 
103: /*=====================================================================*/
104: /*   set external boudary points                                       */
105: /*=====================================================================*/
106: void set_boundary_wall_pot(int n, double u[][MAXDIM], int f[][MAXDIM]){
107: 
108:   int i;
109: 
110:   for(i=0; i<=n; i++){
111:     u[0][i]=0.0;
112:     f[0][i]=1;
113:     
114:     u[n][i]=0.0;
115:     f[n][i]=1;
116: 
117:     u[i][0]=0.0;
118:     f[i][0]=1;
119: 
120:     u[i][n]=0.0;
121:     f[i][n]=1;
122:   }
123:   
124: }
125: 
126: /*=====================================================================*/
127: /*   set circler boundary                                              */
128: /*=====================================================================*/
129: void set_circle(double x0, double y0, double r, double v,
130:         int n, double x[][MAXDIM], double y[][MAXDIM],
131:         double u[][MAXDIM], int f[][MAXDIM]){
132:   
133:   int i, j;
134:   
135:   for(i=0; i<=n; i++){
136:     for(j=0; j<=n; j++){
137:       if((x[i][j]-x0)*(x[i][j]-x0)+(y[i][j]-y0)*(y[i][j]-y0) < r*r){
138:     u[i][j]=v;
139:     f[i][j]=1;
140:       }
141:     }
142:   }
143: }
144: 
145: 
146: /*==========================================================*/
147: /*   3D plot                                                */
148: /*==========================================================*/
149: void plot_3d(FILE *gp, char *data_file,
150:          double rot_x, double rot_z){
151:   
152:   
153:   /* == flowing lines make a graph by using gnuplot == */
154:   
155:   printf("start plot 3D\n");
156:   fprintf(gp, "reset\n");
157:   fprintf(gp, "set terminal postscript eps color\n");
158:   fprintf(gp, "set output \"graph.eps\"\n");
159:   fprintf(gp, "set view %f,%f\n", rot_x, rot_z);
160:   fprintf(gp, "set contour base\n");
161:   fprintf(gp, "set cntrparam levels 20\n");
162:   fprintf(gp, "set hidden3d\n");
163:   fprintf(gp, "set nokey\n");
164:   fprintf(gp, "splot \"%s\" with lines\n", data_file);
165:   fprintf(gp, "set terminal x11\n");
166:   fprintf(gp, "replot\n");
167:   
168: }

---

last update: