Partial differential equations

Heat 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 heat equation  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

p = 0.3;           % coefficient of heat transfer
h = 0.25;          % step in space
a = 0; 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 = xh(2:n-1);    % right hand side (heat sources): f0(x,t) = x
U0 = xh.^2;        % initial condition: u0 = x*x
tau = 0.1;         % step in time
sig = p*tau/(h*h); % coefficient sigma used at finite diff. method
hold off
plot(xh,U0,'g')    % graph of initial condition
hold on
U1(1) = 0; U1(n) = 1;  % boundary conditions - constants

   % computation cycle for time levels
disp('Press a Space tab for the next time level ... ');
for i = 1:9
  U1(2:n-1) = sig*(U0(1:n-2)+U0(3:n))+(1-2*sig)*U0(2:n-1)+tau*f0;
  pause
  plot(xh,U1)
  U0 = U1;
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 = xh.^2;        % initial condition: u0 = x*x
tau = 0.3;         % larger time-step
sig = p*tau/(h*h); % recompute sigma
for i = 1:3
  U1(2:n-1) = sig*(U0(1:n-2)+U0(3:n))+(1-2*sig)*U0(2:n-1)+tau*f0;
  pause;
  plot(xh,U1,'r')
  U0 = U1;
end
disp('The end.');

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

Explanation of the behaviour of the explicit method:

The method can be expressed in matrix form as the Fixed Point Iterations U⁽ⁿ⁺¹⁾ = M U⁽ⁿ⁾ + v⁽ⁿ⁾,
where the iteration matrix M is a Topelitz matrix with values 1-2*sig on the main diagonal and sig on both the first diagonals above and below the main diagonal.

Necessary and sufficient condition for convergence of the FPI is that the spectral radius of M is less than one.

You can check this condition numerically in Matlab for both examples above:

% n = 5; use the value from the previous example
p = 0.3;           % coefficient of heat transfer
h = 0.25;          % spatial step
tau = 0.1;         % time step
sig = p*tau/(h*h); 
I = eye(n); v1 = ones(1,n-1);
M = sig*(diag(v1,1) + diag(v1,-1)) + (1-2*sig)*I;
% or M = toeplitz([1-2*sig , sig, zeros(1,n-2)]);
disp('Spectral radius of iteration matrix is');
disp(max(abs(eig(M)))); % or disp(norm(M));

Heat equation - implicit finite differences

Stability of the implicit method does not depend on the step-length.

Solve the same problem by implicit method and compare the results with the previous ones:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  implicit finite differences for heat equation  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

p = 0.3;           % coefficient of heat transfer
h = 0.25;          % step in space
a = 0; 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 = xh(2:n-1)';   % right hand side (heat sources): f0(x,t) = x (column vector)
U0 = (xh.^2)';     % initial condition: u0 = x*x (column vector)
U1 = U0;           % U1 will contain actual solution
U1(1) = 0; U1(n) = 1;  % boundary conditions - constants
tau = 0.1;         % step in time
sig = p*tau/(h*h); % coefficient sigma used at finite diff. method
hold off
plot(xh,U0,'g')    % graph of initial condition
hold on

   % preparation of the system matrix
A=toeplitz([1+2*sig , -sig, zeros(1,n-4)]);

   % computation cycle for time levels
disp('Press a Space tab for the next time level ... ');
for i = 1:9
    % preparation of the right hand side
  B = U0(2:n-1) + tau*f0;
  B(1) = B(1) + sig * U1(1); B(n-2) = B(n-2) + sig * U1(n); % const. b.c.
    % solving the linear system
  U1(2:n-1) = A\B; 
  pause
  plot(xh,U1)
  U0 = U1;
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 = (xh.^2)';     % initial condition: u0 = x*x (column vector)
tau = 0.3;         % larger time-step
sig = p*tau/(h*h); % recompute sigma

   % recompute the system matrix
A=toeplitz([1+2*sig , -sig, zeros(1,n-4)]);

for i = 1:3
    % preparation of the right hand side
  B = U0(2:n-1) + tau*f0;
  B(1) = B(1) + sig * U1(1); B(n-2) = B(n-2) + sig * U1(n); % const. b.c.
    % solving the linear system
  U1(2:n-1) = A\B; 
  pause;
  plot(xh,U1,'r')
  U0 = U1;
end
disp('The end.');

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