A..5 4_Wave_2D_Contln.c

Wave_2D_Contln.c

   1 /******************************************************************************************
   2  
   3     This program solves a wave equation using the explicit method.
   4  
   5     A wave equation is utt = c * c * uxx.
   6     The function f is the initial shape of wave.
   7  
   8     Time loop in this program consists of 
   9     "Calculation part" and "Draw Part".
  10  
  11     Calculation part is controlled by the interval parameter INTV.
  12     This numerical result is visualized by using "g_contln_2D" of GLSC3D.
  13 
  14 *******************************************************************************************/
  15 
  16 #include<stdio.h>
  17 #include<glsc3d_3.h>
  18 
  19 #define N (100)
  20 #define M (100)
  21 
  22 double f(double x,double y)
  23 {
  24     return 3 * exp(-((x - M_PI / 2)*(x - M_PI / 2) + (y  - 3 * M_PI / 2)*(y  - 3 * M_PI / 2))/ 0.1);
  25 }
  26 
  27 int main()
  28 {
  29     double  x, y;
  30     
  31     double  Lx = 8 * M_PI;
  32     double  Ly = 8 * M_PI;
  33     double  dx = Lx / N;
  34     double  dy = Ly / M;
  35     double  dt = 0.002;
  36     double  c = 1.0;
  37     
  38     double  u0[N + 2][M + 2];
  39     double  u1[N + 2][M + 2];
  40     double  u2[N + 2][M + 2];
  41     
  42     double  lambda_x = c * dt / dx;
  43     double  lambda_y = c * dt / dy;
  44     
  45     int     INTV = 100;
  46     
  47     g_init("Window", 600, 600); //Pixel Size
  48     
  49     g_def_scale_2D(0,                   //ID
  50                    -Lx*0.5, Lx*0.5,     //xmin,xmax
  51                    -Ly*0.5, Ly*0.5,     //ymin,ymax
  52                    20.0, 20.0,          //Window (Left, Top) Position
  53                    560, 560);           //Window Size (x,y)
  54 
  55     //Set initial calculation
  56     {
  57         //Step1
  58         for(int j = 1;j < M + 1; j ++)
  59         {
  60             y = (j - 0.5) * dy;
  61             for(int i = 1;i < N + 1; i ++)
  62             {
  63                 x = (i - 0.5) * dx;
  64                 u0[i][j] = f(x,y);
  65             }
  66         }
  67         //Step2
  68         for(int j = 1;j < M + 1; j ++)
  69         {
  70             u0[0][j] = -u0[1][j];
  71             u0[N + 1][j] = -u0[N][j];
  72         }
  73         for(int i = 0;i < N + 2; i ++)
  74         {
  75             u0[i][0] = -u0[i][1];
  76             u0[i][M + 1] = -u0[i][M];
  77         }
  78         //Step3
  79         for(int j = 1;j < M + 1; j ++)
  80         {
  81             y = (j - 0.5) * dy;
  82             for(int i = 1;i < N + 1; i ++)
  83             {
  84                 x = (i - 0.5) * dx;
  85                 u1[i][j] = u0[i][j] + lambda_x * lambda_x * (u0[i - 1][j] - 2 * u0[i][j] + u0[i + 1][j]) / 2
  86                 + lambda_y * lambda_y * (u0[i][j - 1] - 2 * u0[i][j] + u0[i][j + 1]) / 2;
  87             }
  88         }
  89         //Step4
  90         for(int j = 1;j < M + 1; j ++)
  91         {
  92             u1[0][j] = -u1[1][j];
  93             u1[N + 1][j] = -u1[N][j];
  94         }
  95         for(int i = 0;i < N + 2; i ++)
  96         {
  97             u1[i][0] = -u1[i][1];
  98             u1[i][M + 1] = -u1[i][M];
  99         }
 100     }
 101     
 102     //////////// Start time loop ////////////
 103     for (int i_time = 0; ; i_time++) {
 104 
 105         //////////// Draw Part ////////////
 106         if (i_time%INTV == 0) {
 107             g_cls();                    //Clear window
 108             g_sel_scale(0);             //Select Virtual scale
 109             g_line_color(0.0, 0.0, 0.0, 1.0);
 110             g_boundary();               //Draw Boundary
 111             
 112             g_line_color(1.0, 0.0, 0.0, 1.0);
 113             g_contln_2D(-Lx*0.5 + dx*0.5, Lx*0.5 - dx*0.5, -Ly*0.5 + dy*0.5, Ly*0.5 - dy*0.5, N+2, M+2, u0, 0.01); //Contln_0.01
 114             g_line_color(0.0, 0.0, 1.0, 1.0);
 115             g_contln_2D(-Lx*0.5 + dx*0.5, Lx*0.5 - dx*0.5, -Ly*0.5 + dy*0.5, Ly*0.5 - dy*0.5, N+2, M+2, u0, 0.04); //Contln_0.04
 116             
 117             g_finish();                 //flush Draw buffer
 118             g_sleep(0.01);              //Sleep 0.01 sec
 119         }
 120         
 121         //////////// Calculation Part ////////////
 122         {
 123             //Step5
 124             for(int j = 1;j < M + 1; j ++)
 125             {
 126                 for(int i = 1;i < N + 1; i ++)
 127                 {
 128                     u2[i][j] = 2 * u1[i][j] - u0[i][j]
 129                     + lambda_x * lambda_x *
 130                     (u1[i - 1][j] -2 * u1[i][j] + u1[i + 1][j])
 131                     + lambda_y * lambda_y *
 132                     (u1[i][j - 1] -2 * u1[i][j] + u1[i][j + 1]);
 133                 }
 134             }
 135             //Step6
 136             for(int j = 1;j < M + 1; j ++)
 137             {
 138                 for(int i = 1;i < N + 1; i ++)
 139                 {
 140                     u0[i][j] = u1[i][j];
 141                 }
 142             }
 143             //Step7
 144             for(int j = 1;j < M + 1; j ++)
 145             {
 146                 u0[0][j] = -u0[1][j];
 147                 u0[N + 1][j] = -u0[N][j];
 148             }
 149             for(int i = 0;i < N + 2; i ++)
 150             {
 151                 u0[i][0] = -u0[i][1];
 152                 u0[i][M + 1] = -u0[i][M];
 153             }
 154             //Step8
 155             for(int j = 1;j < M + 1; j ++)
 156             {
 157                 for(int i = 1;i < N + 1; i ++)
 158                 {
 159                     u1[i][j] = u2[i][j];
 160                 }
 161             }
 162             //Step9
 163             for(int j = 1;j < M + 1; j ++)
 164             {
 165                 u1[0][j] = -u1[1][j];
 166                 u1[N + 1][j] = -u1[N][j];
 167             }
 168             for(int i = 0;i < N + 2; i ++)
 169             {
 170                 u1[i][0] = -u1[i][1];
 171                 u1[i][M + 1] = -u1[i][M];
 172             }
 173         }
 174     }
 175     return 0;
 176 }

• 長方形領域上の2次元波動方程式 $\dfrac{1}{c^2}u_{tt}=u_{xx}+u_{yy}$ (境界条件は同次Neumannなのかな) を差分法で解いた解を可視化する (等高線を描く)プログラム。
• g_line_color() も g_area_color() とおなじく、 R, G, B, 不透明度を引数に取る。
• 長方形 $[x_{\text{min}},x_\text{max}]×[y_\text{min},y_\text{max}]$ で定義された関数 $u$ があるとき、定義域を $x$軸方向に $N_x$ 等分, $y$ 軸方向に $N_y$ 等分して格子を作り、格子点での値 $u(x_i,y_j)$ ( $0\le i\le N_x$, $0\le j\le N_y$) を2次元配列 u[nx+1,ny+1] に記録してあるとする。 このとき関数値が level である等高線を描くには、
  g_contln_2D(xmin, xmax, ymin, ymax, nx+1, ny+1, u, level);
  とすれば良い。