1% Leapfrog with the fourth-order Laplacian, otherwise identical to
2% models/leapfrog.m. The only change is `lap4` in place of `lap2`.
3%
4% What it buys is grid dispersion. A discrete Laplacian gets the wave speed
5% slightly wrong at wavelengths that the grid barely resolves, and the error
6% grows with distance travelled: a pulse made of many wavelengths slowly comes
7% apart, its short components lagging behind, leaving a trailing ripple that
8% is entirely numerical. The 5-point stencil's phase error goes as (k*h)^2,
9% the 9-point one's as (k*h)^4, so at the eight or ten points per wavelength
10% where a run is cheap enough to be interactive, the difference is easy to
11% see: raise the frequency until models/leapfrog.m visibly trails, then switch
12% to this one.
13%
14% It is not free. The stencil reads twice as many neighbours, and its
15% stability limit is tighter (c*dt/h <= sqrt(3/8) rather than 1/sqrt(2)), so
16% the app also takes a slightly smaller timestep. Roughly: a little more work
17% per step, in exchange for needing a much coarser grid at the same accuracy.
19function [p, pm, t] = init(npts)
20 p = zeros(npts, 1);
21 pm = zeros(npts, 1);
22 t = zeros(npts, 1);
23end
25function [pn, pold, tn] = step(p, pm, t, x, y, c, sig, dt, f, t0, tw, cw, x0, y0, w, point)
26 u = (t - t0) / tw;
27 env = (1 - cw) * exp(-u .* u) + cw * (0.5 * (1 + tanh(u)));
28 gx = ((x - x0) / w) .^ 2;
29 gy = point * (((y - y0) / w) .^ 2);
30 om = 2*pi*f;
31 s = (om*om) * (env .* sin(om*(t - t0)) .* exp(-(gx + gy)));
33 sd = sig * dt;
34 lap = lap4(p);
35 pn = (2*p - (1 - sd) .* pm + (c*dt).^2 .* lap + (dt*dt) * s) ./ (1 + sd);
37 pold = p;
38 tn = t + dt;
39end