Symmetric Discretization of Linear Boundary Conditions

When solving large-scale PDEs, it's generally preferable (where possible) to have a symmetric discretization so that the conjugate gradient method can be used.  Conjugate gradient uses relatively little memory compared to more general purpose methods (like GMRES) and can also take advantage of symmetric preconditioners (like incomplete LDL^T) that are faster and use lower memory than preconditioners such as incomplete LU.

However applying boundary conditions to the problem can be an issue, since it is easy to end up with a non-symmetric matrix unless the constrained variables are eliminated from the system matrix.  This means that solvers like GMRES have to be used, which is costly in terms of computation and memory.

The following Matlab script gives an example for how to apply linear boundary conditions via an auxiliary boundary condition matrix suggested by a labmate of mine.  It shows both pin (Dirichlet) constraints and gradient (Neumann) boundary conditions at the domain boundary and in the interior.  The only downside is that Neumann constraints are a bit peculiar to specify since one of the variables must be eliminated from the boundary condition matrix.  However this can usually be specified fairly easily, and certainly more easily than eliminating the variable from the system matrix.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Example for showing symmetric discretization  %%
%% of general linear boundary conditions.        %%
%% (c) James Gregson 2013                        %%
%% Please send questions or comments to:         %%
%% james.gregson@gmai.com                        %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% solve a 1D Laplace equaiton with 8 mesh points
% subject to the constraint that the leftmost 
% grid point has value 1.0 and the rightmost 
% grid point has value 2.0 and the fifth is 
% constrained to the value of the fourth minus 0.2
%
% then solves the same problem, but with a biharmonic
% operator, to show that the result is different 
% from just pinning the values; i.e. the biharmonic
% energy is applied to the constrained grid points. 
% If this weren't the case, both problems would have 
% the same solution.

% the approach is to use a standard symmetric 
% discretization of the Laplace operator A for all
% grid points (including those constrained by
% boundary conditions) and a second boundary
% condition matrix that enforces the constraints
% specifically, the boundary conditions are
% represented as:
%
% u = B v + c
%
% where u has the same dimensions as the Laplace
% operator and B has one column for every
% unconstrained grid-point.  c is a vector of
% constants. this allows the problem: 
%
% A u = f
%
% to be written as:
%
% A*B*v = f - A*c
%
% however this is system is rectangular and 
% non-symmetric.  This can be remedied by 
% multiplying through by B', giving:
%
% B'*A*B*v = B'*(f-A*c)
% 
% which is square and symmetric.

% define the objective function, in this case a 
% 6-variable 1D Laplace equation. The Laplacian
% is discretized with a standard 3-point stencil
A = [  1, -1,  0,  0,  0,  0,  0,  0;
      -1,  2, -1,  0,  0,  0,  0,  0;
       0, -1,  2, -1,  0,  0,  0,  0;
       0,  0, -1,  2, -1,  0,  0,  0;
       0,  0,  0, -1,  2, -1,  0,  0;
       0,  0,  0,  0, -1,  2, -1,  0;
       0,  0,  0,  0,  0, -1,  2, -1;
       0,  0,  0,  0,  0,  0, -1,  1 ];

% we are solving the Laplace equation, so the
% right hand side is zero
f = [ 0, 0, 0, 0, 0, 0, 0, 0 ]';

% the boundary condition matrix.  
B = [  0,  0,  0,  0,  0; % first grid point, set to c[1]
       1,  0,  0,  0,  0;
       0,  1,  0,  0,  0;
       0,  0,  1,  0,  0;
       0,  0,  1,  0,  0; % fifth grid point, set to c[4]-0.2
       0,  0,  0,  1,  0;
       0,  0,  0,  0,  1;
       0,  0,  0,  0,  0 ]; % last grid point, set to c[8]

% the boundary condition constant matrix
c = [ 1, 0, 0, 0, -0.2, 0, 0, 2 ]';

% form the system with boundary conditions
M = B'*A*B
rhs = B'*( f - A*c );

% solve for the primary variables
x = M\rhs;

% use the primary variables to solve for the
% grid values
u1 = B * x + c;

% now solve a biharmonic problem with the same constraints
% in the same manner, note the A' factor in the  system
% matrix and right-hand-side.
M = B'*A'*A*B;
rhs = B'*( f - A'*A*c );
x = M\rhs;
u2 = B * x + c;

% plot the results
figure
hold on
plot( u1, '-ro' );
plot( u2, '-b+' );
legend( 'Laplace', 'Biharmonic' );  

The plot produced by this script is shown below:

Both solves reproduce the boundary conditions correctly, but the biharmonic solve finds a solution that is smoother than the Laplace solution around the interior constraints.