Numerical solution to heat equation using FTCS

Hi,

I expect that error should increase as final time T increases but opposite is happening:

Norm of error =  1.1194e-003 at T = 0.9
dt = 0.02

Norm of error =  1.5872e-006 at T = 2
dt = 0.02

Norm of error =  4.9910e-012 at T = 4
dt = 0.02

I have checked my code several times but cannot figure out what is wrong. The code follows:

function [errout, ue] = heatFTCS(dt, nx, k, L, T, errPlots)
% heatFTCS solves 1D heat equation with the FTCS scheme
% Input:  dt - step size in time
%         nx - number of lattice points
%         k - thermal diffusivity
%         L - length of space
%         T - final time
%         errPlots - flag for plotting results
% Output: errout - error in the numerical solution
%         ue - exact solution

% Set default input arguments
if nargin < 1, dt = 0.02; end
if nargin < 2, nx = 11; end
if nargin < 3, k = 1/6; end
if nargin < 4, L = 1; end
if nargin < 5, T = 4.0; end
if nargin < 6, errPlots = 0; end

% determine space grid size, number of time steps and some intermediate values
dx = L/(nx-1);
nt = ceil(T/dt + 1);
r = k*(dt/dx^2);

if r > 0.5
    fprintf('\nr = %g, solution may not converge',r);
end

% create spatial grid and initialize numerical solution
x = linspace(0, L, nx)';
U = zeros(nx, nt);

% Set initial and boundary conditions
U(:,1) = sin(2*pi*x/L);
U(1,:) = 0; U(nx,:) = 0;

% employ numerical scheme
for  m = 2 : nt
    for  i = 2 : nx-1
        U(i,m) = r*U(i-1, m-1) + (1 - 2*r)*U(i, m-1) + r*U(i+1, m-1);
    end
end

% Evaluate exact solution and absolute error
u = sin(2*pi*x/L)*exp(-4*k*(pi/L)^2*T);
err = abs(U(:,nt) - u);

% Specify output arguments, print and plot results
if nargout > 0, errout = err; end
if nargout > 1, ue = u; end

fprintf('\nNorm of error = %12.4e at T = %g\n', norm(err), T);
fprintf('dt = %g\n', dt);

if ~errPlots, return; end
clf reset;
subplot(2,1,1)
plot(x, U(:,nt), '-', x,u, 'o');
legend('FTCS (U)','Exact (u)');
xlabel(sprintf('\\fontname{math}x'),'fontsize',10);
title(sprintf('Numerical (FTCS) vs Exact solution, \\fontname{math}T = %g, dt = %g\n',T,dt),'fontsize',10)
subplot(2,1,2)
plot(x, err,'--r');
xlabel(sprintf('\\fontname{math}x'),'fontsize',10); ylabel(sprintf('\\fontname{math}U - u'),'fontsize',10);
title(sprintf('Error in numerical scheme, \\fontname{math}T = %g, dt = %g\n',T,dt),'fontsize',10)

Thank you!

[Walter Roberson]

Please retag this question after reading the guide to tags, http://www.mathworks.co.uk/matlabcentral/answers/43073-a-guide-to-tags

[Tom]

Why would you expect the error to increase? As far as I can see, you've started with a sinusoidal distribution (symmetric, average being zero), and you're letting it diffuse without any boundary conditions. From this, you'd expect everything to reach a uniform temperature of zero.