波動方程式を解くプログラム

001: #include <stdio.h>
002: #include <unistd.h>
003: 
004: #define NT 800
005: #define NX 200
006: #define SLEEP_TIME 20000
007: 
008: void set_xt (int nx, int nt, double x[], double t[]);
009: void initial_condition(int nx, int nt, double alpha, double u[][NT+1]);
010: /*=====================================================================*/
011: /*   main                                                              */
012: /*=====================================================================*/
013: int main(){
014:   double u[NX+1][NT+1];
015:   double x[NX+1], t[NT+1];
016:   double xmin, xmax, tmin, tmax;
017:   double delta_x, delta_t, alpha;
018:   int i, j;
019:   FILE *gp;
020: 
021:   xmin=0.0;
022:   xmax=1.0;
023: 
024:   tmin=0.0;
025:   tmax=2.7;
026: 
027:   delta_x = (xmax-xmin)/NX;
028:   delta_t = (tmax-tmin)/NT;
029:   alpha = (delta_t/delta_x)*(delta_t/delta_x);
030: 
031:   printf("alpha = %f\n",alpha);
032: 
033:   if(alpha >= 1.0){
034:     printf("alpha must be less than 1. Now alpha = %f\n",alpha);
035:     return 1;
036:   }
037: 
038:   set_xt(NX, NT, x, t);
039:   initial_condition(NX, NT, alpha, u);
040: 
041:   gp = popen("gnuplot -persist","w");
042:   fprintf(gp, "set xrange [-0.1:1.1]\n");
043:   fprintf(gp, "set yrange [-0.6:0.6]\n");
044:   
045:   for(j=0; j<=1; j++){
046:     printf("t=%f\n",t[j]);
047:     fprintf(gp, "plot '-' with lines linetype 1\n");
048:     for(i=0; i<=NX; i++){
049:       fprintf(gp,"%f\t%f\n", x[i], u[i][j]);
050:     }
051:     fprintf(gp, "e\n");
052:     usleep(SLEEP_TIME);
053:   }
054: 
055: 
056:   for(j=2; j<=NT; j++){
057:     printf("%d\tt=%f\n",j,t[j]);
058:     for(i=1; i<NX ; i++){
059:       u[i][j]=2.0*u[i][j-1]-u[i][j-2]
060:     +alpha*(u[i+1][j-1]-2.0*u[i][j-1]+u[i-1][j-1]);
061:     }
062:     fprintf(gp, "plot '-' with lines linetype 1\n");
063:     for(i=0; i<=NX; i++){
064:       fprintf(gp,"%f\t%f\n", x[i], u[i][j]);
065:     }
066:     fprintf(gp,"e\n");
067:     usleep(SLEEP_TIME);
068:   }
069: 
070:   pclose(gp);
071: 
072: 
073:   return 0;
074: }
075: 
076: /*=====================================================================*/
077: /*   set x-coordinate and time                                         */
078: /*=====================================================================*/
079: void
080: set_xt (int nx, int nt, double x[], double t[])
081: {
082: 
083:   int i;
084: 
085:   for (i=0; i<=nx; i++){
086:     x[i]=(double)i/nx;
087:   }
088: 
089:   for (i=0; i<=nt; i++){
090:     t[i]=(double)i/nt;
091:   }
092:   
093: }
094: 
095: 
096: /*=====================================================================*/
097: /*   set initial condition                                             */
098: /*=====================================================================*/
099: void initial_condition(int nx, int nt, double alpha, double u[][NT+1])
100: {
101:   int i;
102: 
103:   /* --------- boundary condition at x=0 and x=1  ------- */
104: 
105:   for(i = 0; i <= nt ; i++){
106:     u[0][i]=0.0;
107:     u[nx][i]=0.0;
108:   }
109: 
110:   /* --------- condition at t=0 ------------ */
111: 
112:   for (i = 1; i < nx ; i++){
113:     if(i <= nx/2){
114:       u[i][0]=(double)i/nx;
115:     }else{
116:       u[i][0]=1.0-(double)i/nx;
117:     }
118:   }
119: 
120:   /* --------- condition at t=delta t ------------ */
121: 
122:   for (i = 1; i < nx ; i++){
123:     u[i][1]=u[i][0]+alpha/2.0*(u[i+1][0]-2*u[i][0]+u[i-1][0]);
124:   }
125: 
126: }

last update: