Partial differential equations

Wave equation - explicit finite differences

Ilustration, how stability of the explicit method depends on the step-length.

Save the following program with an example to a file name.m and run it in Matlab:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  explicit finite differences for wave equation  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

c = 2;             % wave speed
h = 0.2;           % step in space
a = -1; b = 1;     % interval of solution (in space)
xh = a:h:b;        % division of the interval with step h
n = length(xh);    % the number of nodes
f0 = 2*xh(2:n-1);  % right hand side: f0(x,t) = 2*x
U0 = 2*xh.^2;      % initial condition: u0 = 2*x*x
dU0 = (ones(1,n) - xh).*sin(pi*xh/2); % initial condition: du0/dx
hold off
plot(xh,U0,'g')    % graph of initial solution
hold on

tau = 0.08;        % step in time
sig = c*tau/h;     % coefficient sigma used at finite diff. method
s = sig^2;         % all formulas use sigma squared only

   % the first time level computed from initial conditions:
U1(2:n-1) = U0(2:n-1) + tau * dU0(2:n-1);
U1(1) = 2*exp(-tau);  % b.c. at left - time dependent
U1(n) = 2;            % b.c. at right - constant
plot(xh,U1)

   % cycle for time levels
disp('Press a Space tab for the next time level ... ');
for k = 2:8
  U2(2:n-1) = s*(U1(1:n-2)+U1(3:n))+2*(1-s)*U1(2:n-1)-U0(2:n-1)+tau^2*f0;
  U2(1) = 2*exp(-k*tau);  % b.c. at left
  U2(n) = 2;              % b.c. at right
  pause
  plot(xh,U2)
  U0 = U1;
  U1 = U2;
end
disp('The end.');
disp('');

   % change the time step and compute once again
disp('Now the same with larger time-step:'); 
disp('Press a Space tab for the next time level ... ');
U0 = 2*xh.^2;      % initial condition: u0 = 2*x*x
tau = 0.16;        % larger time-step
sig = c*tau/h;     % recompute sigma
s = sig^2;        

   % the first time level computed from initial conditions:
U1(2:n-1) = U0(2:n-1) + tau * dU0(2:n-1);
U1(1) = 2*exp(-tau);  % b.c. at left - time dependent
U1(n) = 2;            % b.c. at right - constant
plot(xh,U1)

for k = 2:4
  U2(2:n-1) = s*(U1(1:n-2)+U1(3:n))+2*(1-s)*U1(2:n-1)-U0(2:n-1)+tau^2*f0;
  U2(1) = 2*exp(-k*tau);  % b.c. at left
  U2(n) = 2;              % b.c. at right
  pause
  plot(xh,U2,'r')
  U0 = U1;
  U1 = U2;
end
disp('The end.');

%%%%%%%%%%%%%%%%%%
%  end of program  
%%%%%%%%%%%%%%%%%%