1% The leapfrog scheme for the wave equation u_tt = c^2 u_xx.
2%
3% Centred second differences in both time and space:
4%
5% u(j, n+1) = 2 u(j, n) - u(j, n-1) + lam^2 (u(j+1, n) - 2 u(j, n) + u(j-1, n)),
6%
7% with lam = c dt / h the Courant number. It is second-order accurate, it is
8% stable exactly when lam <= 1, and it has a property that no other scheme in
9% this tutorial shares: at lam = 1 it is not an approximation at all.
10%
11% u is the displacement at step 0 and uold at step -1, one entry per node.
12% bc is 'fixed', in which case u has N+1 entries and the first and last are
13% held at zero, or 'periodic', in which case it has N entries and node N is
14% followed by node 1. The scheme advances M steps and returns
15%
16% y M-by-1 the displacement at node jr after each step
17% U numel(u)-by-K the whole field after each of the first K steps
18%
19% The inner loop is written out in full, one node at a time; numbl compiles it,
20% and in this form it is about six times faster than the vectorised slice.
22function [y, U] = leapfrog(u, uold, lam, M, bc, jr, K)
23u = u(:);
24uold = uold(:);
25n1 = numel(u);
26l2 = lam ^ 2;
27y = zeros(M, 1);
28U = zeros(n1, K);
29unew = zeros(n1, 1);
30periodic = strcmp(bc, 'periodic');
31for n = 1:M
32 for j = 2:n1 - 1
33 unew(j) = 2 * u(j) - uold(j) + l2 * (u(j + 1) - 2 * u(j) + u(j - 1));
34 end
35 if periodic
36 unew(1) = 2 * u(1) - uold(1) + l2 * (u(2) - 2 * u(1) + u(n1));
37 unew(n1) = 2 * u(n1) - uold(n1) + l2 * (u(1) - 2 * u(n1) + u(n1 - 1));
38 else
39 unew(1) = 0;
40 unew(n1) = 0;
41 end
42 uold = u;
43 u = unew;
44 y(n) = u(jr);
45 if n <= K
46 U(:, n) = u;
47 end
48end
49end