Next: A..7 6_Turing_2D_ColorMap.c Up: A. 講習会のサンプル・プログラムを読む Previous: A..5 4_Wave_2D_Contln.c
| 5_Turing_2D_Contln.c |
1 /******************************************************************************************
2
3 This program solves reacrtion-diffusion equations using the explicit method.
4
5 Reacrtion-diffusion equations are ut = du uxx + f(u,v) and vt = dv vxx + g(u,v).
6
7 Time loop in this program consists of
8 "Calculation part" and "Draw Part".
9
10 Calculation part is controlled by the interval parameter INTV.
11 This numerical result is visualized by using "g_contln_2D" of GLSC3D.
12
13 *******************************************************************************************/
14
15 #include<stdio.h>
16 #include<stdlib.h>
17 #include<glsc3d_3.h>
18
19 #define N (100)
20 #define M (100)
21
22 double f(double u, double v)
23 {
24 return u * (1 - u * u) - v;
25 }
26 double g(double u, double v)
27 {
28 return 3.0 * u - 2.0 * v;
29 // return 3.0 * (u + 0.05) - 2.0 * v;
30 // return 3.0 * (u - 0.05) - 2.0 * v;
31 }
32
33 int main()
34 {
35 double Lx = 100.0;
36 double Ly = 100.0;
37 double dx = Lx / N;
38 double dy = Ly / M;
39 double dt = 0.001;
40 double du = 1.0;
41 double dv = 10.0;
42
43 double lambda_u = du * dt / (dx * dx);
44 double lambda_v = dv * dt / (dy * dy);
45
46 double u0[N + 2][M + 2];
47 double u1[N + 2][M + 2];
48 double v0[N + 2][M + 2];
49 double v1[N + 2][M + 2];
50
51 int INTV = 200;
52
53 g_init("Window", 600, 600); //Pixel Size
54
55 g_def_scale_2D(0, //ID
56 -Lx*0.5, Lx*0.5, //xmin,xmax
57 -Ly*0.5, Ly*0.5, //ymin,ymax
58 20.0, 20.0, //Window (Left, Top) Position
59 560, 560); //Window Size (x,y)
60 //Set initial calculation
61 {
62 //Step1
63 double u_Init_Sum = 0.0;
64 double v_Init_Sum = 0.0;
65
66 unsigned int seed = 2; //Seed for random number generated by rand()
67 srand(seed);
68
69 for(int j = 1;j < M + 1; j ++)
70 {
71 for(int i = 1;i < N + 1; i ++)
72 {
73 double Max = 1.0; //Maximum of random number
74 double Min = 0.0; //Minimum of random number
75 double RandNumber = Min + (int)(rand()*(Max-Min+1.0)/(1.0+RAND_MAX));
76 u0[i][j] = 2.0 * (RandNumber - 0.5);
77 v0[i][j] = 2.0 * (RandNumber - 0.5);
78
79 u_Init_Sum += u0[i][j] * dx * dy;
80 v_Init_Sum += v0[i][j] * dx * dy;
81 }
82 }
83
84 //Step2
85 for(int j = 1;j < M + 1; j ++)
86 {
87 u0[0][j] = u0[N][j];
88 u0[N + 1][j] = u0[1][j];
89 v0[0][j] = v0[N][j];
90 v0[N + 1][j] = v0[1][j];
91 }
92 for(int i = 0;i < N + 2; i ++)
93 {
94 u0[i][0] = u0[i][M];
95 u0[i][M + 1] = u0[i][1];
96 v0[i][0] = v0[i][M];
97 v0[i][M + 1] = v0[i][1];
98 }
99 }
100
101 //////////// Start time loop ////////////
102 for (int i_time = 0; ; i_time++) {
103
104 //////////// Draw Part ////////////
105 if (i_time%INTV == 0) {
106 g_cls(); //Clear window
107 g_sel_scale(0); //Select Virtual scale
108 g_line_color(0.0, 0.0, 0.0, 1.0);
109 g_boundary(); //Draw Boundary
110
111 g_line_color(1.0, 0.0, 0.0, 1.0);
112 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.0); //Contln_0.0
113
114 g_finish(); //flush Draw buffer
115 g_sleep(0.01); //Sleep 0.01 sec
116 }
117
118 //////////// Calculation Part ////////////
119 {
120 //Step3
121 for(int j = 1;j < M + 1; j ++)
122 {
123 for(int i = 1;i < N + 1; i ++)
124 {
125 u1[i][j] = u0[i][j]
126 + lambda_u * (u0[i - 1][j] - 2 * u0[i][j] + u0[i + 1][j])
127 + lambda_u * (u0[i][j - 1] - 2 * u0[i][j] + u0[i][j + 1])
128 + dt * f(u0[i][j], v0[i][j]);
129
130 v1[i][j] = v0[i][j]
131 + lambda_v * (v0[i - 1][j] - 2 * v0[i][j] + v0[i + 1][j])
132 + lambda_v * (v0[i][j - 1] - 2 * v0[i][j] + v0[i][j + 1])
133 + dt * g(u0[i][j], v0[i][j]);
134 }
135 }
136 //Step4
137 for(int j = 1;j < M + 1; j ++)
138 {
139 for(int i = 1;i < N + 1; i ++)
140 {
141 u0[i][j] = u1[i][j];
142 v0[i][j] = v1[i][j];
143 }
144 }
145 //Step5
146 for(int j = 1;j < M + 1; j ++)
147 {
148 u0[0][j] = u0[N][j];
149 u0[N + 1][j] = u0[1][j];
150 v0[0][j] = v0[N][j];
151 v0[N + 1][j] = v0[1][j];
152 }
153 for(int i = 0;i < N + 2; i ++)
154 {
155 u0[i][0] = u0[i][M];
156 u0[i][M + 1] = u0[i][1];
157 v0[i][0] = v0[i][M];
158 v0[i][M + 1] = v0[i][1];
159 }
160 }
161 }
162 return 0;
163 }
|