K. Imposing Dirichlet boundary conditions

There is a simple and widely used method for imposing Dirichlet conditions — direct substitution of the known boundary values into the system (the row-and-column elimination method). Suppose a node ii carries the Dirichlet condition qiq_i =gi= g_i, where gig_i is the known boundary value. Then on the system Aq⃗A \vec{q} =P⃗= \vec{P} we perform the following operations:

  1. Right-hand side modification: For every row jj ≠i\neq i subtract the contribution of the known value qiq_i from the right-hand side: pj:p_j : =pj= p_j −Aji⋅gi- A_{ji} \cdot g_i.
  2. Zeroing the row and column: Replace the ii-th row and ii-th column of the matrix AA with zeros, except for the diagonal element: Aij:A_{ij} : =0,Aji:= 0, \quad A_{ji} : =0,j= 0, \quad j ≠i\neq i.
  3. Setting the diagonal element: Set the diagonal element to one: Aii:A_{ii} : =1= 1.
  4. Setting the boundary value: Replace the ii-th component of the right-hand side with the boundary value: pi:p_i : =gi= g_i.

After these operations the ii-th equation of the system takes the form 1⋅qi1 \cdot q_i =gi= g_i, which guarantees that the boundary condition holds.

Let us return to our one-dimensional example with five nodes q0,q1,q2,q3,q4q_0, q_1, q_2, q_3, q_4. Suppose Dirichlet conditions are prescribed at the boundaries:

q0\displaystyle q_0 =g0,q4\displaystyle = g_0, \quad q_4 =g4.\displaystyle = g_4.

Stationary case

The original system has the form:

a2K⋅q⃗\displaystyle a^2 K \cdot \vec{q} =F⃗.\displaystyle = \vec{F}.

The matrix a2Ka^2 K is full (symmetric and tridiagonal); boundary conditions are imposed by replacing the boundary rows.

Denote AA =a2K= a^2 K and P⃗\vec{P} =F⃗= \vec{F}. We apply the direct method for the nodes q0q_0 and q4q_4.

Step 1: Right-hand side modification for the interior nodes

For the node q1q_1 (index 1):

p1:\displaystyle p_1 : =p1\displaystyle = p_1 −A10⋅g0\displaystyle - A_{10} \cdot g_0 =F1\displaystyle = F_1 −a2⋅(−1l(0)(1))⋅g0\displaystyle - a^2 \cdot \left(-\frac{1}{l_{(0)(1)}}\right) \cdot g_0 =F1\displaystyle = F_1 +a2g0l(0)(1).\displaystyle + \frac{a^2 g_0}{l_{(0)(1)}}.

For the node q3q_3 (index 3):

p3:\displaystyle p_3 : =p3\displaystyle = p_3 −A34⋅g4\displaystyle - A_{34} \cdot g_4 =F3\displaystyle = F_3 −a2⋅(−1l(3)(4))⋅g4\displaystyle - a^2 \cdot \left(-\frac{1}{l_{(3)(4)}}\right) \cdot g_4 =F3\displaystyle = F_3 +a2g4l(3)(4).\displaystyle + \frac{a^2 g_4}{l_{(3)(4)}}.

Steps 2–4: Modifying the matrix and right-hand side for the boundary nodes

For the node q0q_0 (index 0):

  • Zero the first row: A01:A_{01} : =0= 0
  • Zero the first column: A10:A_{10} : =0= 0
  • Set the diagonal: A00:A_{00} : =1= 1
  • Set the right-hand side: p0:p_0 : =g0= g_0

For the node q4q_4 (index 4):

  • Zero the last row: A43:A_{43} : =0= 0
  • Zero the last column: A34:A_{34} : =0= 0
  • Set the diagonal: A44:A_{44} : =1= 1
  • Set the right-hand side: p4:p_4 : =g4= g_4

In the matrix a2Ka^2 K the boundary rows and columns (nodes q0q_0 and q4q_4) are zeroed and the diagonal is set to one; the known values are moved to the right-hand side.

The modified right-hand side vector:

P⃗~\displaystyle \widetilde{\vec{P}} =(g0F1+a2g0l(0)(1)F2F3+a2g4l(3)(4)g4).\displaystyle = \begin{pmatrix} g_0 \\ F_1 + \frac{a^2 g_0}{l_{(0)(1)}} \\ F_2 \\ F_3 + \frac{a^2 g_4}{l_{(3)(4)}} \\ g_4 \end{pmatrix}.

Note: since the source-density function is identically zero on the simplices s1s_1 and s4s_4, we have F0F_0 =0= 0 and F4F_4 =0= 0. However, the components F1,F2,F3F_1, F_2, F_3 may be nonzero.

Non-stationary case

For the non-stationary case the system has the form:

[D+Δt⋅a2K]⋅q⃗n\displaystyle \left[D + \Delta t \cdot a^2 K\right] \cdot \vec{q}_n =Δt⋅F⃗\displaystyle = \Delta t \cdot \vec{F} +D⋅q⃗n−1.\displaystyle + D \cdot \vec{q}_{n-1}.

Denote AA =D= D +Δt⋅a2K+ \Delta t \cdot a^2 K and P⃗\vec{P} =Δt⋅F⃗= \Delta t \cdot \vec{F} +D⋅q⃗n−1+ D \cdot \vec{q}_{n-1}. We apply the same direct method.

Right-hand side modification for the interior nodes

For the node q1q_1:

p1:\displaystyle p_1 : =p1\displaystyle = p_1 −A10⋅g0\displaystyle - A_{10} \cdot g_0 =Δt⋅F1\displaystyle = \Delta t \cdot F_1 +∑jD1jqn−1,j\displaystyle + \sum_j D_{1j} q_{n-1,j} −(l(0)(1)6−Δt⋅a2l(0)(1))⋅g0.\displaystyle - \left(\frac{l_{(0)(1)}}{6} - \frac{\Delta t \cdot a^2}{l_{(0)(1)}}\right) \cdot g_0.

For the node q3q_3:

p3:\displaystyle p_3 : =p3\displaystyle = p_3 −A34⋅g4\displaystyle - A_{34} \cdot g_4 =Δt⋅F3\displaystyle = \Delta t \cdot F_3 +∑jD3jqn−1,j\displaystyle + \sum_j D_{3j} q_{n-1,j} −(l(3)(4)6−Δt⋅a2l(3)(4))⋅g4.\displaystyle - \left(\frac{l_{(3)(4)}}{6} - \frac{\Delta t \cdot a^2}{l_{(3)(4)}}\right) \cdot g_4.

As in the stationary case, the boundary rows and columns of DD +Δt a2K+ \Delta t\, a^2 K are zeroed and the diagonal is set to one, while q0q_0 =g0= g_0 and q4q_4 =g4= g_4 go to the right-hand side.