Boundary conditions

Consider the linear system of equations

A*x = b
where A is a symmetric matrix discretizing a scalar-valued elliptic PDE (A is a discrete Laplace operator). Finding solutions to elliptic PDEs on finite domains requires setting boundary condition e.g., the Dirichlet, Neumann, or periodic.

We consider efficient application of the Dirichlet boundary conditions, i.e., the x value in the boundary nodes is explicitly defined. Assume node with index k is a boundary node. To set the value of x(k)=1 the k‘th row of the system matrix A is modified in the following way

A(k,i)=0 for i=1:n and i!=k
In other words, the diagonal entry of the k‘th row is set to one, and all other entries are set to 0. The k‘th entry of the right-hand side vector b is set to 1

This way the k‘th equation is trivial and reads

Note that after this modification the system matrix A is no longer symmetric. The non-zero columns of the k‘th matrix row have been explicitly set to 0. The symmetry can be regained by setting all the rows in the k‘th matrix column to 0. This is done by eliminating the known x(k) from all the equations it is part of. For every equation i, for which A(i, k) is non-zero the right-hand side vector entry b(k) is modified as

b(k) = b(k) – A(i,k)*x(k)
and the system matrix is modified as

A(i,k) = 0
In general, there are in total n degrees of freedom. Assume there is a set of Dirichlet boundary nodes with indices Bc_ind and boundary values x(Bc_ind) = Bc_val. We denote the indices of the unconstrained nodes as Free_ind. The boundary conditions can be implemented by directly applying the above described modifications to A and b for every boundary node separately.

Free_ind = 1:n;
Free_ind(Bc_ind) = [];
% for all boundary nodes
for id=1:length(Bc_ind)
    % deal with equation k
    k = Bc_ind(id);
    % eliminate the k'th variable
    rows = A(Free_ind,k);
    b(Free_ind) = b(Free_ind) - rows*Bc_val(id);
    % modify A and b
    A(:,k) = 0;
    A(k,:) = 0;
    A(k,k) = 1;
    b(k) = Bc_val(id);

This implementation is very inefficient since every boundary node is processed independently. In particular, reseting row (and column) entries of a sparse matrix to zero introduces costly changes to the sparsity pattern. Avoiding multiple modifications of the sparsity pattern of A is expected to minimize the related overhead.

The code fragment below presents two different approaches to apply the Dirichlet boundary conditions in MATLAB efficiently. Instead of modifying matrix entries, the rows and columns corresponding to the boundary nodes can be all together removed from the system. The update to the right hand side vector is recast into a (dense) matrix-vector multiplication. For symmetric matrices only the upper- or lower-triangular part of A needs to be stored and the implementation is modified accordingly.

 Free = 1:n;
 Free(Bc_ind) = [];
 b = zeros(n,1);
 b(Bc_ind) = Bc_val;
% Approach Ia:
 b = b - (A(:,Bc_ind) + cs_transpose(A(Bc_ind,:)))*Bc_val';
% Approach Ib:
 b = b - (A(:,Bc_ind) + (A(Bc_ind,:))')*Bc_val';
% Approach IIa:
 b = b - A*x - (x'*A)';
% Approach IIb:
 b = b - A*x - A'*x;
% Elimination of trivial equations 1
 A = A(Free_ind,Free_ind);
% Elimination of trivial equations 2
 A(:,Bc_ind) = [];
 A(Bc_ind,:) = [];

Execution time of the above implementations on an example desktop architecture is compared in Figure TODO. The system size is approximately one million degrees of freedom. The number of Bc_ind nodes in the system varies. When the fraction of the boundary nodes is smaller than 1% we observe no substantial difference in the performance of the four methods (Ia,Ib,IIa,IIb) to modify the right-hand side vector b. Most practical two-dimensional models are in this regime. For higher fractions of up to 40% the approaches IIa and IIb perform equally well as for the lower fractions. Implementations Ia and Ib are significantly less efficient in this regime.

Figure TODO presents the execution time of the two studied methods to eliminate the the rows and columns corresponding to the boundary nodes from the system matrix A. When the fraction of the boundary nodes is smaller than 1% the Elimination 1 is slightly more efficient than Elimination 2. For higher fractions of boundary nodes Elimination 1 becomes increasingly more efficient, while the performance of Elimination 2 remains unchanged.

In summary, approaches IIa or IIb should be used to modify the right-hand side vector b, and Elimination 1 should be used to removing the trivial equations for the boundary nodes Bc_ind from the system matrix A.