GMRES: Making the matrix square without solving for boundaries

Question

How do we define the matrix for GMRES, if we do not want to solve the boundary elements but only the interior ones.

I am using pentagonal elements so in a row there are 6 elements (cell itself + 5 fluxes) which means that in a row boundary elements are also accounted. But I need to make the matrix square according to the code here and I do not want to solve for boundary elements? What should I do?

Answer 1

It seems that you are referring to the fluxes that occur at the edge of your domain. Since these are typically known quantities enforced by your boundary conditions, there is no need to solve for those values. There are two common ways to deal with this:

1. Include an equation for each (known) boundary quantity in the form $f_i = F(x_i)$ for $x_i\in\partial\Omega$. This allows for more flexible boundary conditions, but adds an extra equation at each boundary point.
2. Since you know the boundary terms, you may move any quantities involving only these terms to the vector on the right hand side of your matrix equation. This requires a simple change to the right hand side of you equation for each element on the boundary, but cannot be done for some boundary conditions (i.e. Neumann conditions in a finite difference method).

---

As a simple example, consider the following problem:
$$u_{xx} = x^2~~{\rm for}~~x\in[0,1]\\ u(0)=1,~~u(1)=0$$

If this problem is discretized with 2nd order centered finite differences and $\Delta x=0.25$, we get the equations:

\begin{align} u_0 = & 1\\ (u_0 - 2u_{0.25} + u_{0.5} ) / 0.0625 = & 0.0625 \\ (u_{0.25} - 2u_{0.5} + u_{0.75} ) / 0.0625 = & 0.25 \\ (u_{0.5} - 2u_{0.75} + u_{1} ) / 0.0625 = & 0.5625 \\ u_1 = & 0\\ \end{align}

This can either be written in matrix form as 5 equations with 5 unknowns:

$$ \left[\begin{array}{ccccc} 1 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0\\ 0 & 1 & -2 & 1 & 0\\ 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 1 \end{array}\right] \left[\begin{array}{c} u_0\\ u_{0.25}\\ u_{0.5}\\ u_{0.75}\\ u_{1} \end{array}\right] = \left[\begin{array}{c} 1\\ 0.00390625\\ 0.015625\\ 0.03515625\\ 0 \end{array}\right] $$

Or, the known boundary conditions (equations 1 and 5 ) can be used in equations 2 and 4 and the known terms all moved to the right, giving:

$$ \left[\begin{array}{ccc} -2 & 1 & 0\\ 1 & -2 & 1\\ 0 & 1 & -2 \end{array}\right] \left[\begin{array}{c} u_{0.25}\\ u_{0.5}\\ u_{0.75} \end{array}\right] = \left[\begin{array}{c} -0.99609375\\ 0.015625\\ 0.03515625 \end{array}\right] $$

Note:
In the case above, the two systems are exactly equivalent and will give the same answer. However, if I had imposed Neumann conditions instead, the second form would not work since the value at the boundary points would not be known. Instead, the boundary condition needs to be discretized (e.g. by using one-sided differences at the boundary points) and those equations need to be included in the final matrix equation.