Finite difference method with variable coefficient

Question

I want to numerically solve the following elliptic PDE for both when $x \in [0,1]\times [0,1] \subset \mathbb{R}^2$ \begin{align*} -\text{div}(\kappa(x,y)\nabla u(x,y)) &= f(x,y), ~(x,y) \in (0,1)\times(0,1)\\ u(x, 0) &= c, ~\frac{\partial u}{\partial x}(1,y) = 0, ~\frac{\partial u}{\partial x}(0,y) = h(0, y), ~\frac{\partial u}{\partial y}(x,1) = 0, \end{align*} using a finite difference method. When $\kappa$ is constant, this is easy, but I am not sure how to proceed in the case where $\kappa$ is a function. In addition, in this case, how should I incorporate these Neumann boundary conditions?

Answer 1

You can easily apply standard finite differences. Just note that $$ \textrm{div} (k \nabla u) = \nabla k \cdot \nabla u + k \Delta u = k\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial x^2}\right) + \frac{\partial k}{\partial x} \frac{\partial u}{\partial x}+\frac{\partial k}{\partial y} \frac{\partial u}{\partial y}. $$

Your equation on node $(i,j)$ would be

\begin{align*} k_{ij}\left(\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}+\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^2}\right)+\\(k'_x)_{i,j}\frac{u_{i+1,j}-u_{i-1,j}}{2h}+(k'_y)_{i,j}\frac{u_{i,j+1}-u_{i,j-1}}{2h} = -f_{i,j}. \end{align*}