L. Imposing Neumann boundary conditions

Neumann boundary conditions are handled much more simply than Dirichlet ones: the flux prescribed at the boundary immediately gives the vector s⃗\vec{s} (nonzero only at the boundary nodes), which enters only the right-hand side of the system — the stiffness matrix a2Ka^2 K need not be changed.

Stationary case

For the five-node example (the fluxes g0g_0 and g4g_4 are prescribed; the source on the end simplices is zero, F0F_0 =F4= F_4 =0= 0) the modified right-hand side vector is:

P⃗=F⃗+a2s⃗=(a2g0F1F2F3a2g4).\begin{gathered}\vec{P} = \vec{F} + a^2 \vec{s} = \begin{pmatrix} a^2 g_0 \\ F_1 \\ F_2 \\ F_3 \\ a^2 g_4 \end{pmatrix}.\end{gathered}

Non-stationary case

For the non-stationary problem the operator is DD +Δt⋅a2K+ \Delta t \cdot a^2 K, and the flux enters the right-hand side with the factor Δt\Delta t:

[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} +Δt⋅a2s⃗.\displaystyle + \Delta t \cdot a^2 \vec{s}.

The modified right-hand side vector:

P⃗\displaystyle \vec{P} =Δt⋅F⃗\displaystyle = \Delta t \cdot \vec{F} +D⋅q⃗n−1\displaystyle + D \cdot \vec{q}_{n-1} +Δt⋅a2s⃗\displaystyle + \Delta t \cdot a^2 \vec{s} =Δt⋅F⃗\displaystyle = \Delta t \cdot \vec{F} +D⋅q⃗n−1\displaystyle + D \cdot \vec{q}_{n-1} +(Δt⋅a2g0000Δt⋅a2g4).\displaystyle + \begin{pmatrix} \Delta t \cdot a^2 g_0 \\ 0 \\ 0 \\ 0 \\ \Delta t \cdot a^2 g_4 \end{pmatrix}.

The operator DD +Δt⋅a2K+ \Delta t \cdot a^2 K is already nonsingular (the damping matrix DD is positive definite, (D+Δt⋅a2K) 1⃗(D + \Delta t \cdot a^2 K)\,\vec{1} =D1⃗= D \vec{1} ≠0\neq 0), so the non-stationary Neumann problem is solvable even without Dirichlet conditions.