Heat equation via EXPM

Nick Trefethen, October 2010

(Chebfun example pde/Erosion.m)

A well-known PDE problem is the heat equation initial boundary-value problem posed for x in [a,b] and t > 0,

u_t = u_xx,  u(x,0) = u0(x)

with suitable boundary conditions. We can regard this as a time-dependent linear process

u_t = Lu

where L is the operator d^2/dx^2 on [a,b] with the same boundary conditions. The solution is

u(t) = e^(tL) u(0).

In Chebfun we can implement this idea using the EXPM command to compute the operator exponential. Here is an example with Neumann boundary conditions on the interval [0,6]. We start with quite an irregular initial function.

d = [0,6];
u0 = chebfun(@(x) sign((-1).^floor(x.^1.5)),d,'splitting','on');
LW = 'linewidth'; lw = 2; FS = 'fontsize'; fs = 16;
clf, plot(u0,LW,lw), grid on
title(sprintf('t = %4.2f     length = %d',0,length(u0)),FS,fs)
ax = [0 6 -1.2 1.2]; axis(ax)

Here's the solution at t = 0.01. Notice that the narrower spikes have lost more amplitude than the wider ones. The warning message is important: Chebfun does not always give its full accuracy for computations of this kind.

L = chebop(d);          % operator on domain [0,6]
L.op = @(u) diff(u,2);  % 2nd-derivative operator
L.lbc = @(u) diff(u);   % Neumann BC at left
L.rbc = @(u) diff(u);   % Neumann BC at right
dt = 0.01;
expmL = expm(dt*L);     % exponential of the operator
u = expmL*u0;
plot(u,LW,lw), axis(ax), grid on
title(sprintf('t = %4.2f     length = %d',0.01,length(u)),FS,fs)


figure

Warning: Nonsmooth initial data may degrade accuracy in the result. 

Here is the solution at t = 0.02. Now that the function is smooth, there are no further warning messages. The rightmost maximum has extra amplitude, since it effectively corresponded to a wider initial spike thanks to the Neumann boundary condition.

u = expmL*u;
plot(u,LW,lw), axis(ax), grid on
title(sprintf('t = %4.2f     length = %d',0.02,length(u)),FS,fs)

figure

At t = 0.1, there is not much of the original structure left. The length of the chebfun has also been reduced.

u = expm(0.08*L)*u;
plot(u,LW,lw), axis(ax), grid on
title(sprintf('t = %4.2f     length = %d',0.1,length(u)),FS,fs)

Published with MATLAB® 7.14