Method of lines/example implementation
In the main article: The Method of Lines, Part I: Basic Concepts, we discussed some of the basic ideas behind the method of lines (MOL). We now revisit the PDE problem of eqs. (1) to (4) to illustrate the details of constructing a MOL code and to discuss the numerical and graphical output from the code.
The specific initial condition of eq. (2) will be taken as:
| \( u(x,t=0) = sin(\pi x/2)\) | (2) |
The analytical solution to eqs. (1) to (4) ,which we use to evaluate the numerical MOL output, is
| \( u(x,t)=e^{-(\pi^{2}/4)t}sin(\pi x/2) \) | (42) |
Thus, this problem provides us with the means to ascertain precisely how accurate our method is, under well defined circumstances, by permitting the numerical results to be compared directly with the known analytical solution. Verifying the accuracy of the method, when used to solve a representative set of problems, gives confidence that accurate results will be obtained when using the method to solve problems where the solutions are not known apriori.
A main program in Matlab for the MOL solution of eqs. (1) to (4) with the analytical solution, eq. (42), included for comparison with the MOL solution, is listed below.
% File: pde_main.m
% Clear previous files
clear all
clc
%
% Parameters shared with the ODE routine
global ncall ndss
%
% Initial condition
n=21;
for i=1:n
u0(i)=sin((pi/2.0)*(i-1)/(n-1));
end
%
% Independent variable for ODE integration
t0=0.0;
tf=2.5;
tout=linspace(t0,tf,n);
nout=n;
ncall=0;
%
% ODE integration
mf=1;
reltol=1.0e-04; abstol=1.0e-04;
options=odeset('RelTol',reltol,'AbsTol',abstol);
if(mf==1) % explicit FDs
[t,u]=ode15s(@pde_1,tout,u0,options); end
if(mf==2) ndss=4; % ndss = 2, 4, 6, 8 or 10 required
[t,u]=ode15s(@pde_2,tout,u0,options); end
if(mf==3) ndss=44; % ndss = 42, 44, 46, 48 or 50 required
[t,u]=ode15s(@pde_3,tout,u0,options); end
%
% Store numerical and analytical solutions, errors at x = 1/2
n2=(n-1)/2.0+1;
sine=sin(pi/2.0*0.5);
for i=1:nout
u_plot(i)=u(i,n2);
u_anal(i)=exp(-pi^2/4.0*t(i))*sine;
err_plot(i)=u_plot(i)-u_anal(i);
end
%
% Display selected output
fprintf('\n mf = %2d abstol = %8.1e reltol = %8.1e\n',...
mf,abstol,reltol);
fprintf('\n t u(0.5,t) u_anal(0.5,t) err u(0.5,t)\n');
for i=1:5:nout
fprintf('%6.3f%15.6f%15.6f%15.7f\n',...
t(i),u_plot(i),u_anal(i),err_plot(i));
end
fprintf('\n ncall = %4d\n',ncall);
%
% Plot numerical solution and errors at x = 1/2
figure(1);
subplot(1,2,1)
plot(t,u_plot); axis tight
title('u(0.5,t) vs t'); xlabel('t'); ylabel('u(0.5,t)')
subplot(1,2,2)
plot(t,err_plot); axis tight
title('Err u(0.5,t) vs t'); xlabel('t'); ylabel('Err u(0.5,t)')
print -deps pde.eps; print -dps pde.ps
%
% Plot numerical solution in 3D perspective
figure(2);
colormap('Gray');
C=ones(n);
g=linspace(0,1,n); % For distance x
h1 = waterfall(t,g,u',C);
axis('tight');
grid off
xlabel('t, time')
ylabel('x, distance')
zlabel('u(x,t)')
s1 = sprintf('Diffusion Equation - MOL Solution');
sTmp = sprintf('u(x,0) = sin(\\pi x/2 )');
s2 = sprintf('Initial condition: %s', sTmp);
title([{s1}, {s2}], 'fontsize', 12);
v = [0.8616 -0.5076 0.0000 -0.1770
0.3712 0.6301 0.6820 -0.8417
0.3462 0.5876 -0.7313 8.5590
0 0 0 1.0000];
view(v);
rotate3d on;
We can note the following points about this main program:
- After declaring some parameters global so that they can be shared with other routines called via this main program, initial condition (2) is computed over a 21-point grid in \( x \).
% Clear previous files
clear all
clc
%
% Parameters shared with the ODE routine
global ncall ndss
%
% Initial condition
n=21;
for i=1:n
u0(i)=sin((pi/2.0)*(i-1)/(n-1));
end
- The independent variable \( t \) is defined over the interval \( 0 \leq t\leq 2.5 \); again, a 21-point grid is used.
%
% Independent variable for ODE integration
t0=0.0;
tf=2.5;
tout=linspace(t0,tf,n);
nout=n;
ncall=0;
- The 21 ODEs are then integrated by a call to the Matlab integrator ode15s.
mf=1;
reltol=1.0e-04; abstol=1.0e-04;
options=odeset('RelTol',reltol,'AbsTol',abstol);
if(mf==1) % explicit FDs
[t,u]=ode15s(@pde_1,tout,u0,options); end
if(mf==2) ndss=4; % ndss = 2, 4, 6, 8 or 10 required
[t,u]=ode15s(@pde_2,tout,u0,options); end
if(mf==3) ndss=44; % ndss = 42, 44, 46, 48 or 50 required
[t,u]=ode15s(@pde_3,tout,u0,options); end
Three cases are programmed corresponding to mf = 1, 2, 3, for which three different ODEs routines, pde_1 , pde_2 , and pde_3 are called (these routines are discussed subsequently). The variable ndss refers to a library of differentiation routines for use in the MOL solution of PDEs; the use of ndss is illustrated in the subsequent discussion. Note that a stiff integrator, ode15s , was selected because the 21 ODEs are sufficiently stiff that a nonstiff integrator results in a large number of calls to the ODE routine.
- Selected numerical results are stored for subsequent tabular and plotted output.
% Store numerical and analytical solutions, errors at x = 1/2
n2=(n-1)/2.0+1;
sine=sin(pi/2.0*0.5);
for i=1:nout
u_plot(i)=u(i,n2);
u_anal(i)=exp(-pi^2/4.0*t(i))*sine;
err_plot(i)=u_plot(i)-u_anal(i);
end
- Selected tabular numerical output is first displayed.
% Display selected output
fprintf('\n mf = %2d abstol = %8.1e reltol = %8.1e\n',...
mf,abstol,reltol);
fprintf('\n t u(0.5,t) u_anal(0.5,t) err u(0.5,t)\n');
for i=1:5:nout
fprintf('%6.3f%15.6f%15.6f%15.7f\n',...
t(i),u_plot(i),u_anal(i),err_plot(i));
end
fprintf('\n ncall = %4d\n',ncall);
The output from this code is:
mf = 1 abstol = 1.0e-004 reltol = 1.0e-004
t u(0.5,t) u_anal(0.5,t) err u(0.5,t)
0.000 0.707107 0.707107 0.0000000
0.625 0.151387 0.151268 0.0001182
1.250 0.032370 0.032360 0.0000093
1.875 0.006894 0.006923 -0.0000283
2.500 0.001472 0.001481 -0.0000091
ncall = 85
This output indicates that the MOL solution agrees with the analytical solution to at least three significant figures. Also, ode15s calls the derivative routine only 85 times (in contrast with the nonstiff integrator ode45 that requires approximately 5000 - 10000 calls, which clearly indicates the advantage of a stiff integrator for this problem).
- The MOL solution and its error (computed from the analytical solution) are plotted.
% Plot numerical solution and errors at x = 1/2
figure(1);
subplot(1,2,1)
plot(t,u_plot); axis tight
title('u(0.5,t) vs t'); xlabel('t'); ylabel('u(0.5,t)')
subplot(1,2,2)
plot(t,err_plot); axis tight
title('Err u(0.5,t) vs t'); xlabel('t'); ylabel('Err u(0.5,t)')
print -deps pde.eps; print -dps pde.ps
The plotted error output below indicates that the error in the MOL solution varied between approximately \( -3\times 10^{-5} \) and \(16\times 10^{-5}\) which is not quite within the error range specified in the program.
reltol=1.0e-04; abstol=1.0e-04;
The fact that these error tolerances were not satisfied does not necessarily mean that ode15s failed to adjust the integration interval to meet these error tolerances. Rather, the error of approximately \( 1.6 \times 10^{-4} \) is due to the limited accuracy of the second order FD approximation of \( \frac{\partial ^2 u} {\partial x^2} \) programmed in pde_1 . This conclusion is confirmed when the main program calls pde_2 (for mf = 2 ) or pde_3 (for mf = 3 ) as discussed subsequently; these two routines have FD approximations that are more accurate than in pde_1 so the errors fall below the specified tolerances.
This analysis indicates that two sources of errors result from the MOL solution of PDEs such as eq.. (1) :
- errors due to the integration in \( t \) (by ode15s), and
- errors due to the approximation of the spatial derivatives such as \( \frac{\partial ^2 u}{\partial x^2} \) programmed in the derivative routine such as pde_1.
In other words, we have to be attentive to integration errors in the initial and boundary value independent variables.
In summary, a comparison of the numerical and analytical solutions indicates that 21 grid points in \( x \) were not sufficient when using the second order FDs
in pde_1.
However, in general, we will not have an analytical solution such as eq. (42) to determine if the number of spatial grid points is adequate. In this case, some experimentation with the number of grid points, and the observation of the resulting solutions to infer the degree of accuracy or spatial convergence, may be required.
- A 3D plot is also produced.
% Plot numerical solution in 3D perspective
figure(2);
colormap('Gray');
C=ones(n);
g=linspace(0,1,n); % For distance x
h1 = waterfall(t,g,u',C);
axis('tight');
grid off
xlabel('t, time')
ylabel('x, distance')
zlabel('u(x,t)')
s1 = sprintf('Diffusion Equation - MOL Solution');
sTmp = sprintf('u(x,0) = sin(\\pi x/2 )');
s2 = sprintf('Initial condition: %s', sTmp);
title([{s1}, {s2}], 'fontsize', 12);
v = [0.8616 -0.5076 0.0000 -0.1770
0.3712 0.6301 0.6820 -0.8417
0.3462 0.5876 -0.7313 8.5590
0 0 0 1.0000];
view(v);
rotate3d on;
The plotted output below clearly indicates the origin of the lines in the method of lines
The programming of the approximating MOL/ODEs is in one of the three routines called by ode15s. We now consider each of these routines. For mf = 1, pde_1 calls function ut=pde_1(t,u).
% File: pde_1.m
function ut=pde_1(t,u)
%
% Problem parameters
global ncall
xl=0.0;
xu=1.0;
%
% PDE
n=length(u);
dx2=((xu-xl)/(n-1))^2;
for i=1:n
if(i==1) ut(i)=0.0;
elseif(i==n) ut(i)=2.0*(u(i-1)-u(i))/dx2;
else ut(i)=(u(i+1)-2.0*u(i)+u(i-1))/dx2;
end
end
ut=ut';
%
% Increment calls to pde_1
ncall=ncall+1;
We can note the following points about pde_1:
- Some problem parameters are first defined.
function ut=pde_1(t,u)
%
% Problem parameters
global ncall
xl=0.0;
xu=1.0;
xl and xu could have also been set in the main program and passed to pde_1 as global variables. The defining statement at the beginning of pde_1 indicates that the independent variable t and dependent variable vector u are inputs to pde_1, while the output is the vector of t dervatives, ut ; in other words, all of the n ODE derivatives in t must be defined in pde_1.
- The finite difference approximation of eq. (1) is then programmed.
% PDE
n=length(u);
dx2=((xu-xl)/(n-1))^2;
for i=1:n
if(i==1) ut(i)=0.0;
elseif(i==n) ut(i)=2.0*(u(i-1)-u(i))/dx2;
else ut(i)=(u(i+1)-2.0*u(i)+u(i-1))/dx2;
end
end
ut=ut';
The number of ODEs (21) is determined by the length command n=length(u);
so that the programming is general (the number of ODEs can easily be changed in the main program). The square of the FD interval, dx2, is then computed.
- The MOL programming of the 21 ODEs is done in the for loop. For BC (3) , the coding is
if(i==1) ut(i)=0.0;
since the value of \( u(x=0,t) = 0 \) does not change after being set as an initial condition in the main program (and therefore its time derivative is zero).
- For BC (4) , the coding is
elseif(i==n) ut(i)=2.0*(u(i-1)-u(i))/dx2;
which follows directly from the FD approximation of BC (4) ( eq. (37) )
\[ u_x \approx \frac{u(i+1) - u(i-1)}{\Delta x} = 0 \]
or with \( i = n \)
\[ u(n+1) = u(n-1) \]
Note that the fictitious value \( u(n+1) \) can then be replaced in the ODE at \(i=n\) by \(u(n-1) \).
- For the remaining interior points, the programming is
else ut(i)=(u(i+1)-2.0*u(i)+u(i-1))/dx2;
which follows from the FD approximation of the second derivative (eq. (33) )
\[
u_{xx} \approx \frac{(u(i+1) - 2u(i) + u(i-1))}{\Delta x^2}
\]
- Since the Matlab ODE integrators require a column vector of derivatives, a final transpose of ut is required.
ut=ut';
%
% Increment calls to pde_1
ncall=ncall+1;
Finally, the number of calls to pde_1 is incremented so that at the
end of the solution, the value of ncall displayed by the main program
gives an indication of the computational effort required to produce
the entire solution. The numerical and graphical output for this case
(mf=1) was discussed previously.
For mf=2, function pde_2 is called by ode15s.
% File: pde_2.m
function ut=pde_2(t,u)
%
% Problem parameters
global ncall ndss
xl=0.0;
xu=1.0;
%
% BC at x = 0 (Dirichlet)
u(1)=0.0;
%
% Calculate ux
n=length(u);
if (ndss== 2) ux=dss002(xl,xu,n,u); % second order
elseif(ndss== 4) ux=dss004(xl,xu,n,u); % fourth order
elseif(ndss== 6) ux=dss006(xl,xu,n,u); % sixth order
elseif(ndss== 8) ux=dss008(xl,xu,n,u); % eighth order
elseif(ndss==10) ux=dss010(xl,xu,n,u); % tenth order
end
%
% BC at x = 1 (Neumann)
ux(n)=0.0;
%
% Calculate uxx
if (ndss== 2) uxx=dss002(xl,xu,n,ux); % second order
elseif(ndss== 4) uxx=dss004(xl,xu,n,ux); % fourth order
elseif(ndss== 6) uxx=dss006(xl,xu,n,ux); % sixth order
elseif(ndss== 8) uxx=dss008(xl,xu,n,ux); % eighth order
elseif(ndss==10) uxx=dss010(xl,xu,n,ux); % tenth order
end
%
% PDE
ut=uxx';
ut(1)=0.0;
%
% Increment calls to pde_2
ncall=ncall+1;
We can note the following points about pde_2:
- The initial statements are the same as in pde_1. Then the Dirichlet BC at \( x = 0 \) is programmed.
% BC at x = 0 (Dirichlet)
u(1)=0.0;
Acually, the statement u(1)=0.0; has no effect since the dependent variables can only be changed through their derivatives, i.e., ut(1), in the ODE derivative routine. This code was included just to serve as a reminder of the BC at \(x = 0\), which is programmed subsequently.
- The first order spatial derivative \( \frac{\partial u}{\partial x}=u_{x} \), is then computed
% Calculate ux
n=length(u);
if (ndss==2) ux=dss002(xl,xu,n,u); % second order
elseif(ndss== 4) ux=dss004(xl,xu,n,u); % fourth order
elseif(ndss== 6) ux=dss006(xl,xu,n,u); % sixth order
elseif(ndss== 8) ux=dss008(xl,xu,n,u); % eighth order
elseif(ndss==10) ux=dss010(xl,xu,n,u); % tenth order
end
Five library routines, dss002 to dss010, are programmed that use second order to tenth order FD approximations, respectively. Since ndss=4 is specified in the main program, dss004 is used in the calculation of ux.
- BC (4) is then applied followed by the calculation of the second order spatial derivative from the first order spatial derivative.
% BC at x = 1 (Neumann)
ux(n)=0.0;
%
% Calculate uxx
if (ndss== 2) uxx=dss002(xl,xu,n,ux); % second order
elseif(ndss== 4) uxx=dss004(xl,xu,n,ux); % fourth order
elseif(ndss== 6) uxx=dss006(xl,xu,n,ux); % sixth order
elseif(ndss== 8) uxx=dss008(xl,xu,n,ux); % eighth order
elseif(ndss==10) uxx=dss010(xl,xu,n,ux); % tenth order
end
Again, dss004 is called which is the usual procedure (the order of the FD approximation generally is not changed in computing higher order derivatives from lower order derivatives, a process termed stagewise differentiation).
% PDE
ut=uxx';
ut(1)=0.0;
%
% Increment calls to pde_2
ncall=ncall+1;
Note the similarity of the code to the PDE (eq. (1) ), and also the transpose required by ode15s.
The numerical output for this case (mf=2) is:
mf = 2 abstol = 1.0e-004 reltol = 1.0e-004
t u(0.5,t) u_anal(0.5,t) err u(0.5,t)
0.000 0.707107 0.707107 0.0000000
0.625 0.151267 0.151268 -0.0000013
1.250 0.032318 0.032360 -0.0000418
1.875 0.006878 0.006923 -0.0000446
2.500 0.001467 0.001481 -0.0000138
ncall = 62
The plotted error output below indicates that the error in the MOL
solution varied between approximately \(-5\times 10^{-5}\) and
\(3.2\times 10^{-5}\) which is within the error range specified in the
program.
reltol=1.0e-04; abstol=1.0e-04;
Thus, switching from the second order FDs in pde_1 to fourth order
finite differences in pde_2 reduced the spatial truncation error
so that the MOL solution met the specified error tolerances.
For mf = 3, function pde_3 is called by ode15s.
% File: pde_3.m
function ut=pde_3(t,u)
%
% Problem parameters
global ncall ndss
xl=0.0;
xu=1.0;
%
% BC at x = 0
u(1)=0.0;
%
% BC at x = 1
n=length(u);
ux(n)=0.0;
%
% Calculate uxx
nl=1; % Dirichlet
nu=2; % Neumann
if (ndss==42) uxx=dss042(xl,xu,n,u,ux,nl,nu); % second order
elseif(ndss==44) uxx=dss044(xl,xu,n,u,ux,nl,nu); % fourth order
elseif(ndss==46) uxx=dss046(xl,xu,n,u,ux,nl,nu); % sixth order
elseif(ndss==48) uxx=dss048(xl,xu,n,u,ux,nl,nu); % eighth order
elseif(ndss==50) uxx=dss050(xl,xu,n,u,ux,nl,nu); % tenth order
end
%
% PDE
ut=uxx';
ut(1)=0.0;
%
% Increment calls to pde_3
ncall=ncall+1;
We can note the following points about pde_3:
- The initial statements are the same as in pde_1. Then the Dirichlet BC at \(x = 0\) and the Neumann BC at \(x = 1 \) are programmed.
function ut=pde_3(t,u)
%
% Problem parameters
global ncall ndss
xl=0.0;
xu=1.0;
%
% BC at x = 0
u(1)=0.0;
%
% BC at x = 1
n=length(u);
ux(n)=0.0;
Again, the statement u(1)=0.0; has no effect (since the dependent variables can only be changed through their derivatives, i.e., ut(1), in the ODE derivative routine). This code was included just to serve as a reminder of the BC at \( x = 0 \), which is programmed subsequently.
- The second order spatial derivative, \(\frac{\partial^{2}u}{\partial x^{2}}=u_{xx}\), is then computed.
% Calculate uxx
nl=1; % Dirichlet
nu=2; % Neumann
if (ndss==42) uxx=dss042(xl,xu,n,u,ux,nl,nu); % second order
elseif(ndss==44) uxx=dss044(xl,xu,n,u,ux,nl,nu); % fourth order
elseif(ndss==46) uxx=dss046(xl,xu,n,u,ux,nl,nu); % sixth order
elseif(ndss==48) uxx=dss048(xl,xu,n,u,ux,nl,nu); % eighth order
elseif(ndss==50) uxx=dss050(xl,xu,n,u,ux,nl,nu); % tenth order
end
Five library routines, dss042 to dss050, are programmed that use second order to tenth order FD approximations, respectively for a second derivative. Since ndss=44 is specified in the main program, dss0044 is used in the calculation of uxx. Also, these differentiation routines have two parameters that specify the type of BCs:
- nl = 1 or 2 specify a Dirichlet or a Neumann BC, respectively, at the lower boundary value of \( x = xl (= 0)\); in this case, BC (3) is Dirichlet, so nl = 1, and
- nu = 1 or 2 specify a Dirichlet or a Neumann BC, respectively, at the upper boundary value of \(x = xu (= 1)\); in this case, BC (4) is Neumann, so nu = 2.
%
% PDE
ut=uxx';
ut(1)=0.0;
%
% Increment calls to pde_3
ncall=ncall+1;
Again, the transpose is required by ode15s.
The numerical output for this case (mf=3) is:
mf = 3 abstol = 1.0e-004 reltol = 1.0e-004
t u(0.5,t) u_anal(0.5,t) err u(0.5,t)
0.000 0.707107 0.707107 0.0000000
0.625 0.151267 0.151268 -0.0000017
1.250 0.032318 0.032360 -0.0000420
1.875 0.006878 0.006923 -0.0000447
2.500 0.001467 0.001481 -0.0000138
ncall = 62
The plotted error output below indicates that the error in the MOL solution varied between approximately \(-4.8\times 10^{-5}\) and \(3.2\times 10^{-5}\) which is within the error range specified in the program.
reltol=1.0e-04; abstol=1.0e-04;
We conclude this example with the following observation:
As the solution approaches steady state, \(t\rightarrow\infty\),\(u_{t}\rightarrow 0\) and from eq. (1) , \(u_{xx}\rightarrow 0\). As the second derivative vanishes, the solution becomes\[ u_{xx}=0 \]
\( u_{x}=c_{1} \)
\( u=c_{1}x+c_{2} \)
Thus, the steady state solution is linear in \(x\), which can serve as another check on the numerical solution (for BCs (3) and (4) , \(c_{1}=c_{2}=0\) and thus at steady state \(u=0,\) which also follows from the analytical solution, eq. (42) ). This type of special case analysis is often useful in checking a numerical solution. In addition to mathematical conditions such as the linear dependency on \(x\), physical conditions can frequently also be used to check solutions, e.g., conservation of mass, momentum and energy.
Through this example application we have attempted to illustrate the basic steps of MOL/PDE analysis to arrive at a numerical solution of acceptable accuracy. We have also presented some basic ideas for assessing accuracy with respect to time and space (i.e. \( t \) and \( x \)). More advanced applications (e.g., problems expressed as systems of nonlinear PDEs) can be analyzed using the same ideas discussed in this example.
Author William E. Schiesser will be glad to assist by providing codes for other 1D, 2D and 3D PDE problems. In any request please include your name, affiliation and postal mailing address so that we can keep you informed of any changes/additions to the Matlab codes.
Running this example
- Copy and paste the following Matlab routines described above:
- pde_main.m, pde_1.m, pde_2.m and pde_3.m
- Copy and paste the following Matlab library codes:
- Copy all Matlab files into the same directory and run pde_main.m
- To experiment, try changing the various options described above.
