Numerical approximation of the heat equation with Neumann boundary conditions: Method of lines

This article shows how to approximate the heat equation with the method of lines.

1 Heat equation

Heat equation is used to simulate a number of applications related with diffusion processes, as the heat conduction. In the one dimensional case it reads,

[math]u_t(x,t) -u_{xx}(x,t)= f(x,t), \quad x\in(a,b), \quad t\in(0,T) [/math]

[math]u(a,t) = g_1(t), \qquad u_x(b,t)=g_2(t), \quad t\in(0,T) [/math]

[math]u(x,0) = u_0(x), \quad x\in(a,b),[/math]

Here, x is the space variable and t is the time, [math]u(x,t)[/math] represents the temperature of a fine bar located in the segment [math](a,b)[/math], [math]g_1(t)[/math] represent the temperature at the left extreme, [math]g_2(t)[/math] is related with the flux of energy at the right extreme and [math]u_0(x)[/math] the initial distribution of temperatura on the bar at time [math]t=0[/math].

2 Numerical scheme

Space discretization

Take a partition of the space interval [a,b] with step h and denote [math]x_i=a+ih, \; i=0,1,2,...,N[/math], the nodes. Define [math] u_n(t)\sim u(x_n,t),[/math] and replace

[math]-u_{xx}(x,t) \quad \mbox{ by }\quad \frac{-u_{n-1}(t)+2u_n(t) - u_{n+1}(t)}{h^2}[/math]

[math]u_{x}(b,t) \quad \mbox{ by }\quad \frac{u_N+1(t) - u_{N-1}(t)}{2h}[/math]

where [math] u_{N+1}(t),[/math] stands for the solution at the ghost node [math] x_{N+1}=b+h,[/math]. Then, we have the following system of N+2 equations

[math]\left\{ \begin{array}{l} u_n'(t)+\frac{-u_{n-1}(t)+2u_n(t) - u_{n+1}(t)}{h^2}=f(x_n,t), \quad n=1,2,...,{N} \qquad t\gt0, \\ u_0(t)=g_1(t), \qquad t\gt0, \\ {\frac{u_{N+1}(t) - u_{N-1}(t)}{2h}=g_2(t)}, \qquad t\gt0, \\ u_n(0)=u^0(x_n), \qquad n=1,2,...,{N}. \end{array} \right. [/math]

Note that the solution at the node N+1 appears in two equations!

[math] \left\{ \begin{array}{l} {u_0(t)=g_1(t)},\\ u_1'(t)+\frac{{-u_0(t)+}2u_1(t) - u_{2}(t)}{h^2}=f(x_1,t), \\ u_2'(t)+\frac{-u_{1}(t)+2u_2(t) - u_{3}(t)}{h^2}=f(x_2,t), \\ u_3'(t)+\frac{-u_{2}(t)+2u_3(t) - u_{4}(t)}{h^2}=f(x_3,t), \\ ... \\ u_{N}'(t)+\frac{-u_{N-1}(t)+2u_{N}(t){- u_{N+1}(t)}}{h^2}=f(x_N,t), \\ {\frac{u_{N+1}(t)-u_{N-1}(t)}{2h}=g_2(t)}, \end{array} \right. [/math]

Eliminate the equations for [math]u_{0}(t)[/math] and [math]u_{N+1}(t)[/math],

[math] \left\{ \begin{array}{l} u_1'(t)+\frac{2u_1(t) - u_{2}(t)}{h^2}=f(x_1,t)+g_1(t)/h^2, \\ u_2'(t)+\frac{-u_{1}(t)+2u_2(t) - u_{3}(t)}{h^2}=f(x_2,t), \\ u_3'(t)+\frac{-u_{2}(t)+2u_3(t) - u_{4}(t)}{h^2}=f(x_3,t), \\ ... \\ u_{N}'(t)+\frac{-2u_{N-1}(t)+2u_{N}(t)}{h^2}=f(x_N,t)+2g_2(t)/h, \end{array} \right. [/math]

Initial condition:

[math] u_n(0)=u^0(x_n), \qquad n=1,2,...,N. [/math]

In matrix form:

[math] U'(t)+KU(t)=F(t), \qquad {U(0)=U^0} [/math]

where [math] K=\frac1{h^2}\left( \begin{array}{ccccccc} 2& -1& 0& 0& ...& 0 &0 \\ -1& 2& -1& 0& ...& 0 &0 \\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& 0& 0& ...& -2& 2 \end{array} \right), \quad U=\left( \begin{array}{l} u_1\\ u_2 \\... \\u_{N} \end{array} \right) [/math]

[math] U^0=\left( \begin{array}{l} u^0(x_1) \\ u^0(x_2)\\ ... \\ u^0(x_{N}) \end{array} \right), \qquad F=\left( \begin{array}{l} f(x_1,t)+g_1(t)/h^2 \\ f(x_2,t)\\ ... \\ f(x_{N},t)+2g_2(t)/h \end{array} \right) [/math]

Time discretization

We now approximate the solution of the system:

[math] \left\{ \begin{array}{l} U(0)=U^0,\\ U'(t)+KU(t)=F(t), \end{array} \right. [/math]

Chose a time mesh [math]\{t^j \}_{j=0}^J[/math] with time step [math]\Delta t[/math] and define

[math] U^j\sim U(t^j), \qquad f_u(t,U)=-K U+F(t). [/math]

Now we replace the ODE by any numerical method for initial value problems. For example, the Euler method: [math] \left\{ \begin{array}{l} U^0, \\ U^{j+1}= U^j+ \Delta t \; f_u(t^j,U^j) , \qquad j=0,1,...,J-1. \end{array} \right. [/math]

Remark Euler method requires [math]\Delta t \leq h^2/2[/math] for stability.

Other possible methods are:

1. Implicit Euler (unconditionally stable),

[math] \left\{ \begin{array}{l} U^0, \\ U^{j+1}= U^j+ \Delta t \; f_u(t^{j+1},U^{j+1}), \qquad j=0,1,...,J-1 \end{array} \right. [/math]

2. Trapezoidal (or Crank-Nicolson method, unconditionally stable),

[math] \left\{ \begin{array}{l} U^0, \\ U^{j+1}= U^j+ \Delta t \frac{f_u(t^{j},U^{j})+f_u(t^{j+1},U^{j+1})}{2}, \qquad j=0,1,...,J-1 \end{array} \right. [/math]

3. Runge-Kutta

[math] \left\{ \begin{array}{l} U^0, \\ U^{j+1}= U^j+ \frac{\Delta t}{6} \; (k_1+k_2+k_3+k_4), \qquad j=0,1,...,J-1 \\ k_1=f_u(t^{j},U^{j}) \\ k_2=f_u(t^{j}+\frac{\Delta t}{2},U^{j}+\frac{k_1}{2}\Delta t)\\ k_3=f_u(t^{j}+\frac{\Delta t}{2},U^{j}+\frac{k_2}{2} \Delta t)\\ k_4=f_u(t^{j}+\Delta t,U^{j}+k_3 \Delta t) \end{array} \right. [/math]

3 MATLAB code

%% Heat equation approximation: Method of lines-Euler in time
clear all

%% Datos del problema
a=0; b=2; T=1;
% Dato inicial y condiciones frontera
f=@(x,t) (x+t);
u0=@(x) (exp(2*x)-1).*(2-x);
g1=@(t) sin(t);
g2=@(t) sin(t);

%% Discretización
N=20;   % number of space subintervals
mu=0.5; % relation between space and time discretization
h=(b-a)/N; % space step
xi=a+h:h:b; % interior space nodes (only interior nodes have equations with Dirichlet BC)
x=a:h:b;  % space nodes (including boundary ones)
dt=mu*h^2; % time step
t=0:dt:T;  % vector of times

%% Escribimos el sistema: U'=-K*U+F(t), U(0)=U0. Calculamos K,F(t) y U0:
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=K/h^2; % aproximation of -u_xx: matrix N x N
G1=zeros(N,1); G1(1)=1; % to include the boundary term in the first equation
G2=zeros(N,1); G2(N)=1; % % to include the boundary term in the last equation 
F=@(t) f(xi,t)'+g1(t)/h^2*G1+2*g2(t)/h*G2; % vector F(t)
U0=u0(xi)'; % vector with the initial condition

%% Métodos Explícitos: calculamos la función fu para escribir el sistema: U'=fu(t,U)
fu=@(t,U) -K*U+F(t); 

%% Bucle en tiempo
sol(1,:)=u0(x); % Vamos guardando la solución en la matriz sol (N+1)x(length(t))
U=U0; % Creamos un vector U que va a tener la aproximación en cada instante de tiempo. Inicialmente U=U0.
for k=1:length(t)-1
    U=U+dt*fu(t(k),U); % Método de Euler en tiempo 
    % Pare resolver por el método del Trapecio descomentar la línea siguiente:
    % U=(eye(N)+dt*K/2)\((eye(N)-dt*K/2)*U+dt*(F(t(k))+F(t(k+1)))/2);
    sol(k+1,:)=[g1(t(k+1)),U']; % guardamos la solución en la fila correspondiente de la matriz sol.
end

%% Dibujamos
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)

4 Results

Consider the particular case: [math](a,b)=(0,2),\; u_0(x)=(e^{2x}-1)(2-x), \; g_1(t)=g_2(t)=\sin(t), \; f(x,t)=x+t, h=1/10 [/math]

left

--Carlos Castro (discusión) 15:09 31 ene 2013 (CET)