1: import graph3;
2: import palette;
3: import animate;
4:
5: settings.tex = "pdflatex";
6: settings.render = 0;
7: settings.prc = false;
8:
9: int n = 26;
10: real L = 2.5;
11: real dx = 2 * L / n;
12: real CFL = 0.125;
13: real dt = CFL * dx ^ 2;
14:
15: real[][] Z = new real[n][n];
16: real[][] W = new real[n][n];
17:
18: void advanceZ(int iter = 20)
19: {
20: for(int k = 0; k < iter; ++k)
21: {
22: for (int i = 0; i < n; ++i)
23: {
24: for (int j = 0; j < n; ++j)
25: {
26: W[i][j] = 0;
27: }
28: }
29:
30: for (int i = 1; i < n - 1; ++i)
31: {
32: for (int j = 1; j < n - 1; ++j)
33: {
34: W[i][j] = Z[i][j] + (dt / dx ^ 2) *
35: (Z[i+1][j]+Z[i-1][j]+Z[i][j-1]+Z[i][j+1]-4*Z[i][j]);
36: }
37: }
38:
39: for (int i = 0; i < n; ++i)
40: {
41: for (int j = 0; j < n; ++j)
42: {
43: Z[i][j] = W[i][j];
44: }
45: }
46: }
47: }
48:
49: unitsize(1cm);
50: pen[] Pal = Rainbow(96);
51: animation a;
52: surface sf;
53: int endframe = 40;
54:
55: currentprojection = perspective(-20, -18, 18);
56: currentlight = light(1, 1, 10);
57:
58: guide initcond1 = shift((-1, -1))* scale(0.5)*unitcircle;
59: guide initcond2 = shift((0.5, 0))*yscale(1.2)*unitcircle;
60:
61: real f(pair p)
62: {
63: return (inside(initcond1,p) || inside(initcond2,p)) ? 2 : 0;
64: }
65:
66: real g(pair t)
67: {
68: int i = round((n / 2) * (t.x / L + 1));
69: int j = round((n / 2) * (t.y / L + 1));
70:
71: if (i > n - 1) { i = n - 1; }
72: if (j > n - 1) { j = n - 1; }
73:
74: return Z[i][j];
75: }
76:
77: // initialize
78: for( int i = 0; i < n; ++i )
79: {
80: for (int j = 0; j < n; ++j)
81: {
82: Z[i][j] = f((-L, -L) + (2 * L / n) * (i,j));
83: }
84: }
85:
86: for(int fr = 0; fr < endframe; ++fr)
87: {
88: if(fr == 0) // smoothing of initial state; no Spline, but full grid
89: {
90: advanceZ(3);
91: sf = surface(g, (-L, -L), (L, L), nx = n);
92: }
93: else // use Spline and fewer grid points to save memory
94: {
95: sf = surface(g, (-L, -L), (L, L), nx = round(n/2), Spline);
96: }
97:
98: sf.colors(palette(sf.map(zpart), Pal[0:round(48 * max(sf).z)]));
99:
100: draw(sf);
101: a.add();
102: erase();
103: advanceZ(30);
104: }
105:
106: label(a.pdf (delay = 400, "controls,loop"));
107: shipout(bbox(3mm, darkblue + 3bp + miterjoin, FillDraw(fillpen=paleblue)), "pdf");