function [  ] = leapfrog(  )

% The leap-frog scheme is unstable for lambda = 1
lambda = 0.8;
h = 1/40;
M = round(4/h);                         % round to the nearest integer
k = lambda*h;                           %
N = round(6/k);                       % round to the nearest integer
u0 = setInitialCondition(M);            % Set initial condition
u = LeapFrog (M, N, u0, lambda);        % Apply the finite difference scheme
% draw animation
figure(1); clf;
F(N) = struct('cdata',[],'colormap',[]);
x = linspace(-1,3,M+1);
for l = 1:N
    plot(x, u(l,:));
    axis([0,3, -0.2, 1]);
    drawnow;
    F(l) = getframe;
end
movie(F);

end


function [u0] = setInitialCondition(M)
% Set an (M+1)-dim row vector containing the initial condition
x = linspace(-1,3,M+1);        % Define (M+1)-dim row vector containing x coordinates
u0 = zeros(1,M+1);             % Allocate memory for row vector u0
ind = (x>=-0.5) & (x <= 0.5);  % Indices of grid points in [-0.5, 0.5]
u0(ind) = cos(pi*x(ind)).^2;   % Set non-zero values for u0
end

function [u] = LeapFrog(M, N, u0, lambda)
% Implementation of the leap frog scheme
u = zeros(N+1, M+1);
u(1,:) = u0;
% Use forward-time central space to compute the solution at t_1
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);
% The leap frog scheme starts from t_2
for n=3:(N+1)
    u(n, 2:M) = u(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);
end
end

Published with MATLAB® R2014b