d8098edAn interactive tutorial on the wave equation in acoustics: strings, drums and rooms, measured by earJeremy Magland 1% Leapfrog in a rigid rectangular room: the seven-point Laplacian in three
2% dimensions, stable for lam <= 1/sqrt(3).
3%
4% The grid has n(1) by n(2) by n(3) nodes, with the walls on the outermost
5% nodes. A rigid wall means du/dn = 0, which is imposed by mirroring the
6% first node inside each wall into a ghost node just outside it. That
7% mirror is the discrete version of the method of images: the grid's
8% solution in the room is exactly the grid's free-space solution for the
9% source and all its mirror images, so the only error left is the scheme's
10% own dispersion.
11%
12% The room starts silent. At every step the source node receives src(n),
13% already scaled by the driver to dt^2 f(t_n) / h^3, which is a point source
14% of strength f(t) spread over one cell. The result is the pressure at the
15% receiver node after each step.
16%
17% The field is stored as one long vector and the neighbours of node q are
18% q +- 1, q +- sy and q +- sz. numbl runs this form about seven times faster
19% than the same loop written with three subscripts.
21function y = room3d(n, lam, M, is, ir, src)
22nx = n(1) + 2;
23ny = n(2) + 2;
24nz = n(3) + 2;
25sy = nx;
26sz = nx * ny;
27u = zeros(nx * ny * nz, 1);
28up = u;
29un = u;
30l2 = lam ^ 2;
31qs = is(1) + 1 + is(2) * sy + is(3) * sz; % padded, so every index is one further in
32qr = ir(1) + 1 + ir(2) * sy + ir(3) * sz;
33y = zeros(M, 1);
34for step = 1:M
35 % ghosts: the mirror image of the first interior node, on all six faces
36 for c = 1:nz
37 for j = 1:ny
38 b = (c - 1) * sz + (j - 1) * sy;
39 u(b + 1) = u(b + 3);
40 u(b + nx) = u(b + nx - 2);
41 end
42 end
43 for c = 1:nz
44 for i = 1:nx
45 b = (c - 1) * sz + i;
46 u(b) = u(b + 2 * sy);
47 u(b + (ny - 1) * sy) = u(b + (ny - 3) * sy);
48 end
49 end
50 for j = 1:ny
51 for i = 1:nx
52 b = (j - 1) * sy + i;
53 u(b) = u(b + 2 * sz);
54 u(b + (nz - 1) * sz) = u(b + (nz - 3) * sz);
55 end
56 end
57 for c = 2:nz - 1
58 for j = 2:ny - 1
59 b = (c - 1) * sz + (j - 1) * sy;
60 for i = 2:nx - 1
61 q = b + i;
62 un(q) = 2 * u(q) - up(q) + l2 * (u(q + 1) + u(q - 1) + u(q + sy) + u(q - sy) ...
63 + u(q + sz) + u(q - sz) - 6 * u(q));
64 end
65 end
66 end
67 un(qs) = un(qs) + src(step);
68 up = u;
69 u = un;
70 y(step) = u(qr);
71end
72end