function [ ] = leapfrogRAfilter( )
lambda = 0.8;
h = 1/20;
M = round(4/h);
k = lambda*h;
N = round(2.4/k);
u0 = setInitialCondition(M);
[u1, w1] = LeapFrogRA (M, N, u0, lambda, 0);
[u2, w2] = LeapFrogRA (M, N, u0, lambda, 0.1);
[u3, w3] = LeapFrogRA (M, N, u0, lambda, 0.2);
x = linspace(-1,3,M+1);
plot(x, w1(N-1,:), x, w2(N-1,:), x, w3(N-1,:));
legend('\nu = 0', '\nu = 0.1', '\nu = 0.2');
end
function [u0] = setInitialCondition(M)
x = linspace(-1,3,M+1);
u0 = zeros(1,M+1);
ind = (x>=-0.5) & (x <= 0.5);
u0(ind) = cos(pi*x(ind)).^2;
end
function [u, w] = LeapFrogRA(M, N, u0, lambda, nu)
u = zeros(N+1, M+1);
u(1,:) = u0;
w = u;
u(2, 2:M) = u(1,2:M) - 0.5*lambda*u(1,3:(M+1)) + 0.5*lambda*u(1,1:(M-1));
u(2,1) = 0;
u(2,M+1) = u(2,M);
for n=3:(N+1)
u(n, 2:M) = w(n-2,2:M) - lambda*( u(n-1,3:(M+1)) - u(n-1,1:(M-1)) );
u(n,1) = 0;
u(n,M+1) = u(n,M);
w(n-1,:) = u(n-1,:) + nu/2*( u(n,:) - 2*u(n-1,:) + w(n-2,:) );
end
end