d8098edAn interactive tutorial on the wave equation in acoustics: strings, drums and rooms, measured by earJeremy Magland 1% Three explicit schemes for the square membrane.
2%
3% 'fd2' leapfrog in time, the five-point Laplacian in space: second order
4% in both, stable for lam <= 1/sqrt(2). Chapter 4's scheme.
5% 'fd4' leapfrog in time, the fourth-order Laplacian in space, whose 1-D
6% stencil is (-1, 16, -30, 16, -1)/12: fourth order in space but
7% still second order in time, and stable only for lam <= 0.612.
8% 'fd44' fourth order in both, by the modified-equation trick (Dablain
9% 1986): the leapfrog update's leading time error is
10% (dt^2/12) u_tttt = (dt^2/12) c^4 Lap^2 u, so add it back, using
11% the five-point Laplacian applied twice for Lap^2, since the term
12% is already multiplied by dt^2. One more stencil pass per step,
13% and the time step is 1/sqrt(2) again.
14%
15% The wide stencils need values one node beyond the rim. With u = 0 on the
16% rim, the wave equation makes every even derivative of u vanish there too,
17% so the solution continues across the rim as an odd reflection, and that is
18% exactly what the ghost nodes are filled with. For the sine modes of the
19% square the reflection is not an approximation at all.
20%
21% The contract is Chapter 4's: u and uold are (N+1)-by-(N+1), rim included;
22% the result is the displacement at node (ir, jr) after each of M steps, and
23% the whole field after each of the first K.
25function [y, U] = drum2d(u, uold, lam, M, scheme, ir, jr, K)
26n1 = size(u, 1);
27g = n1 + 4; % two ghost layers on every side
28P = zeros(g, g);
29Q = zeros(g, g);
30R = zeros(g, g);
31W = zeros(g, g);
32P(3:n1 + 2, 3:n1 + 2) = u;
33Q(3:n1 + 2, 3:n1 + 2) = uold;
34y = zeros(M, 1);
35U = zeros(n1 * n1, K);
36if strcmp(scheme, 'fd2')
37 s = 0;
38elseif strcmp(scheme, 'fd4')
39 s = 1;
40else
41 s = 2;
42end
43l2 = lam ^ 2;
44l4 = lam ^ 4 / 12;
45a = n1 + 2; % the index of the far rim
46for n = 1:M
47 % odd reflection across the rim, two layers deep
48 for j = 1:g
49 P(2, j) = -P(4, j);
50 P(1, j) = -P(5, j);
51 P(a + 1, j) = -P(a - 1, j);
52 P(a + 2, j) = -P(a - 2, j);
53 end
54 for i = 1:g
55 P(i, 2) = -P(i, 4);
56 P(i, 1) = -P(i, 5);
57 P(i, a + 1) = -P(i, a - 1);
58 P(i, a + 2) = -P(i, a - 2);
59 end
60 if s == 0
61 for j = 4:a - 1
62 for i = 4:a - 1
63 R(i, j) = 2 * P(i, j) - Q(i, j) ...
64 + l2 * (P(i + 1, j) + P(i - 1, j) + P(i, j + 1) + P(i, j - 1) - 4 * P(i, j));
65 end
66 end
67 else
68 if s == 2
69 % W = five-point Laplacian of P (times h^2), zero on the rim,
70 % reflected oddly one layer out
71 for j = 4:a - 1
72 for i = 4:a - 1
73 W(i, j) = P(i + 1, j) + P(i - 1, j) + P(i, j + 1) + P(i, j - 1) - 4 * P(i, j);
74 end
75 end
76 for j = 1:g
77 W(2, j) = -W(4, j);
78 W(a + 1, j) = -W(a - 1, j);
79 end
80 for i = 1:g
81 W(i, 2) = -W(i, 4);
82 W(i, a + 1) = -W(i, a - 1);
83 end
84 end
85 for j = 4:a - 1
86 for i = 4:a - 1
87 L4 = (-P(i + 2, j) + 16 * P(i + 1, j) - 30 * P(i, j) + 16 * P(i - 1, j) - P(i - 2, j) ...
88 - P(i, j + 2) + 16 * P(i, j + 1) - 30 * P(i, j) + 16 * P(i, j - 1) - P(i, j - 2)) / 12;
89 v = 2 * P(i, j) - Q(i, j) + l2 * L4;
90 if s == 2
91 v = v + l4 * (W(i + 1, j) + W(i - 1, j) + W(i, j + 1) + W(i, j - 1) - 4 * W(i, j));
92 end
93 R(i, j) = v;
94 end
95 end
96 end
97 Q = P;
98 P = R;
99 y(n) = P(ir + 2, jr + 2);
100 if n <= K
101 U(:, n) = reshape(P(3:a, 3:a), n1 * n1, 1);
102 end
103end
104end