第 3 章 他の領域のプログラム
3.3 Cassini の oval
Cassiniのovalとは、方程式
((x−x0)2+ (y−y0)2+a2)2 =b4+ 4a2(x−x0)2
で与えられる曲線です。この曲線は、2点(x0+a, y0),(x0−a, y0)からの距離の積が一定値b2 に等しい点の軌跡です。
x, yについて解くと、それぞれ、
x=x0±
√
a2−(y−y0)2±√
b4−4a2(y−y0)2, y=y0±√
−a2−(x−x0)2±√
b4 + 4a2(x−x0)2 となります。ただし、√
内は正でなければならないので、
y=y0±
√
−a2−(x−x0)2+√
b4+ 4a2(x−x0)2
となります。またCassiniのovalはa, bの大小関係で以下のように形を変えます。
a > bのとき (aとbの値が大きいほど左寄り)
a=bのとき (8の字型)
a < bのとき (aとbの値が大きいほど右寄り)
プログラムは以下のものになります。
1 /*
2 * sotsuronC.c --- Cassininのovalでの熱方程式を SW 近似で解く
3 */
45 #include <stdio.h>
6 #include <math.h>
78 /* to use matrix, new_matrix() */
9 #include <matrix.h>
1011 /* to use GLSC */
12 #define G_DOUBLE 13 #include <glsc.h>
1415
16 /* 領域の内部・境界・外部の判定用 */
17 double F(double x, double y, double a, double b, double x0, double y0){
18 return b*b
19 - sqrt(((x-x0)*(x-x0)+(y-y0)*(y-y0)+a*a)*((x-x0)*(x-x0)+(y-y0)*(y-y0)+a*a) 20 -4.0*a*a*(x-x0)*(x-x0));
21 }
2223 /* 東側の境界 */
24 double E(double x, double y, double a, double b, double x0, double y0){
25 if(x <= x0){
26 return x0-sqrt(a*a - (y-y0)*(y-y0) - sqrt(b*b*b*b-4*a*a*(y-y0)*(y-y0)));
27 }
28 else return x0+sqrt(a*a - (y-y0)*(y-y0) + sqrt(b*b*b*b-4*a*a*(y-y0)*(y-y0)));
29 } 30
31 /* 西側の境界 */
32 double W(double x, double y, double a, double b, double x0, double y0){
33 if(x <= x0){
34 return x0-sqrt(a*a - (y-y0)*(y-y0) + sqrt(b*b*b*b-4*a*a*(y-y0)*(y-y0)));
35 }
36 else return x0+sqrt(a*a - (y-y0)*(y-y0) - sqrt(b*b*b*b-4*a*a*(y-y0)*(y-y0)));
37 }
3839 /* 北側の境界 */
40 double N(double x, double a, double b, double x0, double y0){
41 return y0+sqrt(-a*a - (x-x0)*(x-x0) + sqrt(b*b*b*b+4*a*a*(x-x0)*(x-x0)));
42 }
4344 /* 南側の境界 */
45 double S(double x, double a, double b, double x0, double y0){
46 return y0-sqrt(-a*a - (x-x0)*(x-x0) + sqrt(b*b*b*b+4*a*a*(x-x0)*(x-x0)));
47 } 48
49 /* 初期条件 */
50 double f(double x, double y, double a, double b, double x0, double y0){
51 return b*b
52 - sqrt(((x-x0)*(x-x0)+(y-y0)*(y-y0)+a*a)*((x-x0)*(x-x0) 53 +(y-y0)*(y-y0)+a*a) -4.0*a*a*(x-x0)*(x-x0));
54 }
55 /* ここまで */
5657 double min(double x, double y){
58 if(x < y){
59 return x;
60 }
61 else return y;
62 } 63
6465 int main() 66 {
67 int Nx, Ny, i, j, n, skip, nMax;
68 double hx, hy, lambda, tau, Tmax, t, dt, ew, es, ee, en, er, mine, maxtau;
69 double x, y, hw, he, hn, hs;
70 double xi, xiw, xie, yj, yjs, yjn;
71 double uw,ue,un,us;
72 double a, b, A, xmin, xmax, ymin, ymax, x0, y0, distance;
73 matrix u, newu;
7475
76 /* 領域のx軸、y軸の長さ */
77 printf("xmin= "); scanf("%lf",&xmin);
78 printf("xmax= "); scanf("%lf",&xmax);
79 printf("ymin= "); scanf("%lf",&ymin);
80 printf("ymax= "); scanf("%lf",&ymax);
8182 /* 分割数 */
83 printf("Nx= "); scanf("%d",&Nx);
84 printf("Ny= "); scanf("%d",&Ny);
8586 /* 円の中心 */
87 printf("x0= "); scanf("%lf",&x0);
88 printf("y0= "); scanf("%lf",&y0);
8990 /* 円の半径 */
91 printf("a= "); scanf("%lf",&a);
92 printf("b= "); scanf("%lf",&b);
93
94 /* 倍率 */
95 printf("A= "); scanf("%lf",&A);
96
97 /* 距離 */
98 printf("distance= "); scanf("%lf",&distance);
99
100 /* 格子間の長さ */
101 hx = (xmax-xmin)/Nx;
102 hy = (ymax-ymin)/Ny;
103 printf("hx=%g, hy=%g\n", hx, hy);
104105 /* 誤差 */
106 er= 1.0e-14;
107108 mine= 1.0;
109 maxtau= 1.0;
110111 /****************最小εと最大τを求める *****************/
112113 for (i = 0;i <= Nx; i++){
114 xi = xmin+i*hx;
115 xiw = xmin+(i-1)*hx;
116 xie = xmin+(i+1)*hx;
117
118 for (j = 0;j <= Ny; j++){
119 yj = ymin+j*hy;
120 yjs = ymin+(j-1)*hy;
121 yjn = ymin+(j+1)*hy;
122123 if(F(xi, yj, a, b, x0, y0) >= er){
124 /* WEST */
125 if(F(xiw, yj, a, b, x0, y0) >= er){
126 hw= hx;
127 ew= 1.0;
128 }
129 else if(fabs(F(xiw, yj, a, b, x0, y0)) < er){
130 ew= 1.0;
131 hw= hx;
132 }
133 else{
134 hw= xi - W(xi, yj, a, b, x0, y0);
135 ew= hw/hx;
136 } 137 /* SOUTH */
138 if(F(xi, yjs, a, b, x0, y0) >= er){
139 hs= hy;
140 es= 1.0;
141 }
142 else if(fabs(F(xi, yjs, a, b, x0, y0)) < er){
143 hs= hy;
144 es= 1.0;
145 }
146 else{
147 hs= yj - S(xi, a, b, x0, y0);
148 es= hs/hy;
149 }
150 /* EAST */
151 if(F(xie, yj, a, b, x0, y0) >= er){
152 he= hx;
153 ee= 1.0;
154 }
155 else if(fabs(F(xie, yj, a, b, x0, y0)) < er){
156 ee= 1.0;
157 he= hx;
158 }
159 else{
160 he= E(xi, yj, a, b, x0, y0) - xi;
161 ee= he/hx;
162 }
163 /* NORTH */
164 if(F(xi, yjn, a, b, x0, y0) >= er){
165 en= 1.0;
166 hn= hy;
167 }
168 else if(fabs(F(xi, yjn, a, b, x0, y0)) < er){
169 en= 1.0;
170 hn= hy;
171 }
172 else{
173 hn= N(xi, a, b, x0, y0) - yj;
174 en= hn/hy;
175 }
176 }
177 else {
178 ew = ee = en = es = 1.0;
179 hw = he = hx;
180 hn = hs = hy;
181 }
182183 mine= min(mine, min(min(ew,ee),min(en,es)));
184 maxtau= min(maxtau,
185 ew*es*ee*en*hx*hx*hy*hy/(2.0*(hx*hx*ew*ee+hy*hy*es*en)));
186
187 }
188 }
189
190 printf("εの最小値= %g\n", mine);
191 printf("τの最大値=%g\n", maxtau);
192 if(mine < er){
193 printf("mine=%g: *******εが小さすぎます********\n", mine);
194 }
195
196 if ((u = new_matrix(Nx + 1, Ny + 1)) == NULL) { 197 fprintf(stderr, "配列 u を確保できませんでした。");
198 exit(1);
199 }
200 if ((newu = new_matrix(Nx + 1, Ny + 1)) == NULL) { 201 fprintf(stderr, "配列 newu を確保できませんでした。");
202 exit(1);
203 }
204
205 printf("Tmax= "); scanf("%lf", &Tmax);
206207 printf("τ (≦ %g )== ", maxtau); scanf("%lf", &tau);
208209 lambda = tau/(hx*hx) + tau/(hy*hy);
210 printf("λ= %g になりました。\n", lambda);
211212 printf("Δt= "); scanf("%lf", &dt);
213
214 skip = rint(dt / tau);
215 if (skip == 0) {
216 printf("Δ tが小さすぎるので、Δt= τ とします。\n");
217 skip = 1;
218 dt = skip * tau;
219 }
220221 g_init("Meta", 250.0, 170.0);
222 g_device(G_BOTH);
223 g_def_scale(0, xmin, xmax, ymin, ymax, 100.0, 100.0, 100.0, 100.0);
224 g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
225 g_sel_scale(0);
226
227 for (i = 0; i <= Nx; i++) 228 for (j = 0; j <= Ny; j++)
229 u[i][j] = f(xi, yj, a, b, x0, y0)/A;
230231 nMax = rint(Tmax / tau);
232
233 for (n = 0; n <= nMax; n++) {
234235 /**********************領域***********************/
236237 for (i = 0; i <= Nx; i++){
238 xi = xmin+i*hx;
239 xiw = xmin+(i-1)*hx;
240 xie = xmin+(i+1)*hx;
241242 for (j = 0; j <= Ny; j++){
243 yj = ymin+j*hy;
244 yjn = ymin+(j+1)*hy;
245 yjs = ymin+(j-1)*hy;
246
247 if(F(xi, yj, a, b, x0, y0) >= er){
248 /* WEST */
249 if(F(xiw, yj, a, b, x0, y0) >= er){
250 hw= hx;
251 uw= u[i-1][j];
252 }
253 else if(fabs(F(xiw, yj, a, b, x0, y0)) < er){
254 hw= hx;
255 uw= u[i-1][j];
256 }
257 else{
258 hw= xi - W(xi, yj, a, b, x0, y0);
259 uw= 0;
260 }
261 /* EAST */
262 if(F(xie, yj, a, b, x0, y0) >= er){
263 he=hx;
264 ue=u[i+1][j];
265 }
266 else if(fabs(F(xie, yj, a, b, x0, y0)) < er){
267 he=hx;
268 ue=u[i+1][j];
269 }
270 else{
271 he= E(xi, yj, a, b, x0, y0) - xi;
272 ue= 0;
273 }
274 /* NORTH */
275 if(F(xi, yjn, a, b, x0, y0) >= er){
276 hn=hy;
277 un=u[i][j+1];
278 }
279 else if(fabs(F(xi, yjn, a, b, x0, y0)) < er){
280 hn=hy;
281 un=u[i][j+1];
282 }
283 else{
284 hn= N(xi, a, b, x0, y0) - yj;
285 un= 0;
286 }
287 /* SOUTH */
288 if(F(xi, yjs, a, b, x0, y0) >= er){
289 hs=hy;
290 us=u[i][j-1];
291 }
292 else if(fabs(F(xi, yjs, a, b, x0, y0)) < er){
293 hs=hy;
294 us=u[i][j-1];
295 }
296 else{
297 hs= yj - S(xi, a, b, x0, y0);
298 us= 0;
299 }
300 }
301 else{
302 u[i][j]= 0;
303 hw= he= hx;
304 hn= hs= hy;
305 uw= ue= un= us= 0;
306 }
307308 newu[i][j] = u[i][j]
309 + 2.0*tau*((ue-u[i][j])/he - (u[i][j]-uw)/hw) /(he+hw) 310 + 2.0*tau*((un-u[i][j])/hn - (u[i][j]-us)/hs )/(hn+hs);
311312 }
313 }
314315 for (i = 0; i <= Nx; i++) 316 for (j = 0; j <= Ny; j++) 317 u[i][j] = newu[i][j];
318 if (n % skip == 0){
319 g_cls();
320 g_hidden2(xmax-xmin, ymax-ymin, 10.0, 0.0, -10.0, distance, 321 10.0, 20.0, 50.0, 50.0, 150.0, 100.0, u, Nx+1, Ny+1,
322 1, G_SIDE_NONE, 1, 1);
323 }
324 t = tau * n;
325 }
326 /* マウスでクリックされるのを待つ */
327 g_sleep(-1.0);
328 /* ウィンドウを消す */
329 g_term();
330 return 0;
331 }