How and when to impose inhomogeneous Dirichlet boundary conditions in Newton method solver for a PDE?

Question

I am applying a central Finite Difference scheme in space and an implicit Euler scheme in time on a variant of the 2d Burgers equation, of the form:$$u_t + uu_x + uu_y = \nu(u_{xx} + u_{yy})$$ where $u$ is an unknown variable that depends on x,y and t, and $u_t,u_x,u_y$ are the single derivatives of $u$ with respect to t,x and y respectively, and $u_{xx},u_{yy}$ the double derivatives of $u$ with respect to x and y. For this problem, I have inhomogeneous Dirichlet boundary conditions that vary with time. After discretisation, I apply the Newton method to the $u$ vector(calculating Jacobian matrix and value of function at every node for right hand side) at every time step and solve a linear system of equations. The thing is, I do not understand whether I should enforce boundary conditions or not.

My current approach is to enforce boundary conditions at the end of each Newton iteration, that is after I update $u$, and then at the end of each time step. This feels a bit like I am "cheating".

This is a sentiment I noticed in Method of lines for inhomogeneous Dirichlet conditions, but there the approach involves Finite Element method. Sadly, I have very little knowledge in this method and I could not translate the answer there to my problem properly. So to summarise: Is my approach correct? If not, what should usually be done?

Answer 1

For Dirichlet BCs, this process is probably a bit simpler than you are thinking. I will demonstrate the procedure for a 1D Poisson problem, but it translates very easily to other problems with structured grids in any dimension.

Consider the problem $$-u'' =f \ \ \forall x\in (-0,1), \quad u(0)=a, \ u(1)=b$$ with the standard 3-point finite difference stencil. Denote the discretization of $u$ and $f$ as $U$ and $F$, respectively. The finite difference equations then read $$U_{i-1} - 2U_i + U_{i+1} = h^2F_i, \ \ i=1\dots N.$$ Notice here that I am indexing such that $U_1$ through $U_N$ correspond to interior points of the domain, where the PDE is explicitly defined. The boundary conditions are given by $U_0 = a$ and $U_{N+1} = b$. As written, this is a system of $N+2$ equations and $N+2$ unknowns, but these extra two "unknowns" will be eliminated soon.

Now consider the finite difference equation at the node adjacent to the boundary, i.e., $i=1$ and $i=N$. There, we have $$U_0 - 2U_1 + U_2 = h^2F_1 \ \ \text{and} \ \ U_{N-1} - 2U_N + U_{N+1} = h^2F_N.$$ However, notice that we already know $U_0$ and $U_{N+1}$ exactly, so we can substitute and rearrange these equations to obtain $$- 2U_1 + U_2 = h^2F_1 - a \ \ \text{and} \ \ U_{N-1} - 2U_N = h^2F_N - b,$$ respectively. We now have $N$ equations with $N$ unknowns given by the system $AU=\tilde{F}$, where $A\in\mathbb{R}^{N\times N}$ is tridiagonal with diagonals $1/h^2$, $-2/h^2$, and $1/h^2$, $U\in\mathbb{R}^N$, and $\tilde{F}\in\mathbb{R}^N$ is an augemented load vector given by $$\tilde{F} = F + \begin{pmatrix}-\frac{a}{h^2} \\0 \\ \vdots \\ 0 \\\frac{-b}{h^2}\end{pmatrix}.$$

The key point here is that you do some work on paper to completely eliminate the "unkowns" corresonding to the boundary nodes, then the boundary conditions become part of the residual for the linear/Newton solve. This has the effect of ensuring consistency of your boundary conditions throughout the solution, i.e., your interior nodes are actually evolving according to the true BC and they aren't using the wrong BC internally then being set to the right BC, but the interior nodes are wrong. This also often has the effect of making the Jacobian of the problem more scructured since now the only unkowns being considered are ones subject to the same dynamics, being the PDE. For nonlinear problems, once typically isn't able to completely isolate everything corresponding to the boundary condition on the RHS as in this problem, but doing so should make your code more robust and easier to generalize to other Dirichlet BCs.

Answer 2

(I am not a mathematician, so there might be some inaccuracies)

To bring your question in a broader context, note the following:

Your system can be written as $$L(u,u',u'',x,t)=0$$ on $\Omega$ with in-homogeneous linear boundary operator $B$ $$B(u,u',x,t)=g(x,t)$$ on $\partial\Omega$.

Now we make use of the linearity of the boundary operator and define the solution $u=v+w$ as a superposition of two different solutions. Be $w$ the solution of the in-homogeneous boundary problem only $$B(w,w',x,t)=g(x,t)$$ and $v$ the solution of the homogeneous problem $B(v,v',x,t)=0$ and a modified in-homogeneous system $\hat{L}$

$$\begin{eqnarray}L(u,u',...,x,t)&=&0\\ L(v+w,(v+w)',...,x,t)&=&0\\ \hat{L}(v,v',...,x,t) &=& f(w,w',...,x,t) \end{eqnarray}$$

Now, how does this transfer to your Newton method?

Consider the implementation of the Newton iteration with homogeneous boundary conditions only, but with an additional right hand side of your PDE.

Then solve in 3 steps:

1. Solve for $w$
2. Construct $\hat{L}$ with $f$
3. Solve with Newton iteration for $v$

return $u=v+w$.

Since the extension into the domain $\Omega$ of $w$ does not really matter as long as you can easily construct $\hat{L}$ and $f$, you can e.g. choose $w=0$ for all interior nodes.

My current approach is to enforce boundary conditions at the end of each Newton iteration, that is after I update u, and then at the end of each time step. This feels a bit like I am "cheating".

Now with the slight change of your iteration technique, the cheating is actually mathematically sound and is nothing else but $u=v+w$ and $w=0$ on interior nodes.

Actually, I believe there are probably some compatibility conditions in the choice of function spaces for $w$ and $v$ since there need to exist a trace operator for $w$ on $\Omega$. Maybe a mathematician can expand on that.