Numerical Inversion of Elliptic Operator

Question

I am trying to solve the following elliptic differential equation.

$$ \frac{\partial^2\psi}{\partial R^2} + \frac{\partial^2\psi}{\partial Z^2} - \frac{1}{R}\frac{\partial\psi}{\partial R} = S(R,Z) $$

The interior points are discretized in a N-1 x N-1 grid as

$$ \left(A_{RR}+A_{ZZ}-\frac{A_R}{R_{i,j}}\right)\psi_{i,j}=S_{i,j}-b^{RR}_{i,j}-b^{ZZ}_{i,j} + \frac{b^R_{i,j}}{R_{i,j}} $$

where

$$ A_{RR}=I\otimes A_{RR}^{1D} $$ $$ A_{RR}^{1D} = \frac{1}{\Delta R^2}\left[\begin{array}{cccc} -2 & 1\\ 1 & -2 & 1\\ & 1 & -2 & 1\\ & & 1 & -2 \end{array}\right] $$

$$ A_{ZZ}=A_{RR}^{1D} \otimes I $$

$$ A_{R}=I \otimes A_R^{1D} $$

$$ A_R^{1D}=\frac{1}{2\Delta R}\left[\begin{array}{cccc} 0 & 1\\ -1 & 0 & 1\\ & -1 & 0 & 1\\ & & -1 & 0 \end{array}\right] $$

where I am assuming a 4x4 interior grid for simplicity, and the spacing is equal to $ \Delta R$ in both directions. The other terms on the right hand side are vectors containing the boundary conditions. I have written the following code:

def invertOpDirichletB(S, R, bcMat):
    dR = R[1] - R[0]
    N = len(R) - 2
    ROr = R
    R = R[1:-1]

    # central difference matrix, R direction
    eR = np.ones(N)
    CD_data = np.array([-eR, eR]) / R
    CD_diags = np.array([-1, 1])
    CD = sparse.diags(CD_data, CD_diags, shape=(N, N))
    CD2D = sparse.kron(sparse.eye(N), CD)

    # laplacian matrix
    DRR_data = np.array([eR, -2*eR, eR])
    DRR_diags = np.array([-1, 0, 1])
    DRR = sparse.diags(DRR_data, DRR_diags, shape=(N, N))

    eZ = np.ones(N)
    DZZ_data = np.array([eZ, -2*eZ, eZ])
    DZZ_diags = np.array([-1, 0, 1])
    DZZ = sparse.diags(DZZ_data, DZZ_diags, shape=(N, N))

    lap = sparse.kron(sparse.eye(N), DRR) + sparse.kron(DZZ, sparse.eye(N))

    LHS = lap / dR ** 2 - 0.5 * CD2D / dR

    Sint = S[1:-1, 1:-1]

    bcRMat = np.zeros([N, N])
    bcRMat[:, -1] = (bcMat[:, -1])[1:-1]
    bcRMat[:, 0] = (bcMat[:, 0])[1:-1]
    bcRRMat = bcRMat / dR ** 2
    bcRMat[:, 0] = 0.5 * bcRMat[:, 0] / (ROr[0] * dR)
    bcRMat[:, -1] = 0.5 * bcRMat[:, -1] / (ROr[-1] * dR)


    bcZZMat = np.zeros([N, N])
    bcZZMat[0] = (bcMat[0])[1:-1]
    bcZZMat[-1] = (bcMat[-1])[1:-1]
    bcZZMat /= dR ** 2
    bcMatT = -bcRRMat - bcZZMat + bcRMat

    RHS = np.ravel(Sint) + np.ravel(bcMatT)

    soln = sparse.linalg.spsolve(LHS, RHS)
    
    soln = soln.reshape([N, N])
    soln = np.pad(soln, 1, constant_values = 0)
    soln[0] = bcMat[0]
    soln[-1] = bcMat[-1]
    soln[:, -1] = bcMat[:, -1]
    soln[:, 0] = bcMat[:, 0]
    
    return soln

Which works fine with homogeneous boundary conditions, for a source term like this one. When I set the boundary conditions to 1, however, the result is not quite right. Here's what happens when I apply the operator again, which should return the source term. There is clearly something wrong with my boundary condition on the left but I cannot figure out what it is.

Any help would be very much appreciated.

Thank you

Answer 1

In case it's useful for anyone later, I fixed my problem by changing the line bcRMat[:, 0] = 0.5 * bcRMat[:, 0] / (ROr[0] * dR) to bcRMat[:, 0] = -0.5 * bcRMat[:, 0] / (ROr[0] * dR).