Approximate Neumann BVP operator by a matrix

Question

I am considering approximating the operator of the BVP \begin{cases} -u''+u'=g,&\quad x\in [a,b]\\ u'(a)=-1, u'(b)=1, \end{cases} by a matrix.

---

I tried to use the idea of finite difference method. I first consider the operator of the homogenuous BVP \begin{cases} -u''+u'=\lambda u,&\quad x\in [a,b]\\ u'(a)=u'(b)=0. \end{cases} Suppose $[a,b]$ is divided into uniform subintervals and each subinterval $(x_i,x_{i+1})$is represented by the midpoint $\bar{x}_i$, mesh size $h=x_{i+1}-x_{i}$. The first order differential operator is approximated by $$\dfrac{\partial }{\partial x}|_{\frac{1}{2}}=\dfrac{1}{2}(\dfrac{\partial }{\partial x}|_{0}+\dfrac{\partial }{\partial x}|_{1})=\dfrac{1}{2}(0+\dfrac{1 }{h}(u_{\frac{3}{2}}-u_{\frac{1}{2}})),$$ where we use the Neumann BC. Similarly, we get $$\dfrac{\partial }{\partial x}|_{N-\frac{1}{2}}=\dfrac{1}{2h}(-u_{N-\frac{3}{2}}+u_{N-\frac{1}{2}}),$$ for the last element. In this way, we can approximate the $\frac{\partial}{\partial x}$ of the BVP by the matrix $$D=\begin{bmatrix} -1 & 1 & & & \\ -1 & 0 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 0 & 1 \\ & & & -1 & 1 \end{bmatrix}.$$

For the second order operator, $$\dfrac{\partial^2 }{\partial x^2}|_{\frac{1}{2}}=\dfrac{1}{h}(\dfrac{\partial }{\partial x}|_{1}-\dfrac{\partial }{\partial x}|_{0})=\dfrac{1}{h^2}(u_{\frac{3}{2}}-u_{\frac{1}{2}}),$$ for the first element. In this way, we can approximate the $\frac{\partial^2}{\partial x^2}$ of the BVP by the matrix $$DD=\begin{bmatrix} -1 & 1 & & & \\ 1 & -2 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -2 & 1 \\ & & & 1 & -1 \end{bmatrix}.$$

But for the inhomogenous case, we have $$\dfrac{\partial }{\partial x}|_{\frac{1}{2}}=\dfrac{1}{2}(\dfrac{\partial }{\partial x}|_{0}+\dfrac{\partial }{\partial x}|_{1})=\dfrac{1}{2}(-1+\dfrac{1 }{h}(u_{\frac{3}{2}}-u_{\frac{1}{2}})),$$ then I don't know how to handle the constant, since it's independent of $u$. Is there anyway to include the BC into the matrix?

Answer 1

If your domain is defined as $a = x_0 < x_1 < \dotsb < x_N = b$ you can add ghost points $$ \begin{aligned}x_{-1} = a-\Delta x\,,\\ x_{N+1} = b+\Delta x\,,\end{aligned} $$ then impose the boundary conditions as $$ \begin{aligned} u'(a) = \frac{u_1 - u_{-1}}{2\Delta x}\,, \\ u'(b) = \frac{u_{N+1} - u_{N-1}}{2\Delta x}\,. \end{aligned} $$ Rearranging for $u_{-1}$ and $u_{N+1}$, you can express the finite difference operator by changing the matrix and including a boundary correction vector.

For example, $$ \begin{aligned} u''(a) &= \frac{u_1 -2u_0 + u_{-1}}{\Delta x^2}\\ &= \frac{u_1 - 2u_0}{\Delta x^2} + \frac{u_1 - 2\Delta x u'(a)}{\Delta x^2}\\ &= \frac{2u_1 - 2u_0}{\Delta x^2} - \frac{2u'(a)}{\Delta x}\,. \end{aligned}$$

Answer 2

An answer to your question "how to handle the constant" is that it has to be moved to the right hand side of your linear algebraic system. In another words this part of inhomogeneous Neumann boundary condition shall NOT be included to the matrix, because, as you mentioned, it is independent of $u$. It will cause the only non-zero term in the right hand side of your system of linear algebraic equations if you discretize your homogeneous BVP.