Numerical approximation of the heat equation with Dirichlet 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(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),g_2(t)[/math] represent the temperature at the extremes 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 [0,L] with step h. 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]
Then, we have the system
[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-1 \qquad t\gt0, \\ u_0(t)=g_1(t), \qquad t\gt0, \\ u_{N}(t)=g_2(t), \qquad t\gt0, \\ u_n(0)=u^0(x_n), \qquad n=1,2,...,N-1. \end{array} \right. [/math]
The first N+1 equations can be written as
[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-1}'(t)+\frac{-u_{N-2}(t)+2u_{N-1}(t) {- u_{N}(t)}}{h^2}=f(x_{N-1},t), \\ {u_N(t)=g_2(t)}, \end{array} \right. [/math]
Eliminate the equations for [math]u_0(t)[/math] and [math]u_{N}(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-1}'(t)+\frac{-u_{N-2}(t)+2u_{N-1}(t)}{h^2}=f(x_{N-1},t)+g_2(t)/h^2, \end{array} \right. [/math]
The initial condition is tranlated into the N-1 equations:
[math] u_n(0)=u^0(x_n), \qquad n=1,2,...,N-1. [/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& ...& -1& 2 \end{array} \right), \quad U=\left( \begin{array}{l} u_1\\ u_2 \\... \\u_{N-1} \end{array} \right) [/math]
[math] U^0=\left( \begin{array}{l} u^0(x_1) \\ u^0(x_2)\\ ... \\ u^0(x_{N-1}) \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-1},t)+g_2(t)/h^2 \end{array} \right) [/math]
Time discretization
[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-h; % 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, U0 y F
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=K/h^2; % aproximation of -u_xx: matrix (N-1) x (N-1)
U0=u0(xi)'; % vector with the initial condition
G1=zeros(N-1,1); G1(1)=1; % to include the boundary term in the first equation
G2=zeros(N-1,1); G2(N-1)=1; % % to include the boundary term in the last equation
F=@(t) f(xi,t)'+g1(t)/h^2*G1+g2(t)/h^2*G2; % vector F
%% 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-1)+dt*K/2)\((eye(N-1)-dt*K/2)*U+dt*(F(t(k))+F(t(k+1)))/2);
sol(k+1,:)=[g1(t(k+1)),U',g2(t(k+1))]; % 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]
--Carlos Castro (discusión) 15:09 31 ene 2013 (CET)
