Is the following PDE boundary value problem well-posed? My Question

*

*Is the following Poisson boundary value problem well-posed, as stated?

*

*If so, how could I go about solving it?

*If not, what would it need to be well-posed?



*Does it satisfy the "compatibility constraint" mentioned in the Update section below?

*(Bonus!) Are there any potential issues for solving this numerically (e.g. using an implicit finite difference scheme to solve for $\Phi$)?

My Approaches Thus Far
I have tried solving this numerically in MATLAB and Mathematica. Both programs claim to have solved the problem, but upon inspection, neither the governing equation nor any of the boundary conditions are enforced!
Here's the problem:
It's an axisymmetric Poisson problem in spherical coordinates.
Domain

$0 < R_{min} \le r \le R_{max}$
$0 < \theta_{min} \le \theta \le \theta_{max} < \pi$

Governing Equation

$2 r \frac{\partial\Phi}{\partial r} + r^2 \frac{\partial^2\Phi}{\partial r^2} + \cot \theta \frac{\partial\Phi}{\partial\theta} + \frac{\partial\Phi}{\partial\theta^2}  = \frac{A}{8 \pi r^2} \csc\theta \left(1+3 \cos(2 \theta)\right) $

The constant $A$ is known.
Radial Boundary Conditions

$ \frac{\partial\Phi}{\partial r} \bigg|_{r \rightarrow R_{min}}  = -\frac{A \sin \theta}{4 \pi R_{min}^3}$
$ \frac{\partial\Phi}{\partial r} \bigg|_{r \rightarrow R_{max}}  = -\frac{A \sin \theta}{4 \pi R_{max}^3}$

Latitudinal Boundary Conditions

$\frac{\partial\Phi}{\partial\theta} \bigg|_{\theta \rightarrow \theta_{min}}  = \frac{A \cos \left(\theta_{min}\right)}{2 \pi r^2}$
$\frac{\partial\Phi}{\partial\theta} \bigg|_{\theta \rightarrow \theta_{max}}  = \frac{A \cos \left(\theta_{max}\right)}{2 \pi r^2}$

Additional Constraint(?)
Since $\Phi(r,\theta)$ is only constrained by its derivatives, its solution is only uniquely determined up to an additive constant. That is, since $\vec{\nabla} \Phi = \vec{\nabla} \left(\Phi + C\right)$, where $C$ is a uniform constant, we need to specify a value for $\Phi(r,\theta)$ at some point we'll call ($r^\star, \theta^\star$):

$\Phi(r^\star, \theta^\star) = C$

in order to get a unique solution. The particular value of $C$ is irrelevant for me, since I'm only interested in the gradient of $\Phi$ anyway.
Background

*

*The equations given above are greatly-simplified versions of a more general problem I am tackling in my research.

*Physically, I'm using the electric scalar potential to calculate the electric field in the given domain, which is an 2D (axisymmetric) slice of a sphere.

*Later, it would be neat for the domain to extend to the poles, but for simplicity (and for numerics' sake!), I want the domain to stop short of the poles.

Update!
Evidently, the problem should satisfy a compatibility constraint. See §1.3 of this document for details. I'm not sure how it applies to this PDE, but I would like this to be addressed in the answer.
Thanks!
 A: Your equation is
$$
     L\Phi=\frac{\partial}{\partial r}r^{2}\frac{\partial\Phi}{\partial r}
    + \frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\sin\theta\frac{\partial\Phi}{\partial\theta}=\frac{A}{8\pi r^{2}}\frac{1}{\sin\theta}(1+3\cos(2\theta)).
$$
Try to find a solution $\Phi_{h}$ of the above equation with homogenous boundary conditions, and then solve $L\Phi_{j}=0$ with homogenous boundary conditions, except for one of the $j=1,2,3,4$ conditions. That would give a solution $\Phi=\Phi_{h}+\Phi_{1}+\Phi_{2}+\Phi_{3}+\Phi_{4}$. The solutions $\Phi_j$ should be unique, and separation of variables should work for these if you're careful about the weight functions you're using.
A: Short answer
It should be well-posed for any domain satisfying $0<R_{min} < R_{max}$ and $0 \le \theta_{min} < \theta_{max} \le \pi$. 
Unfortunately, and probably due to some error, the answer seems to depend on whether I analyze this problem in 3D (assuming axisymmetry) or in 2D (on a plane of constant $\phi$)....
Long answer
To be well-posed, the problem must satisfy the consistency condition (or "compatibility constraint"), which is derived thus:
Derivation of the Compatibility Constraint for a general Poisson equation with purely Neumann boundary conditions:
Given the PDE and its boundary conditions:
\begin{align}
    \nabla^2 \Phi \;=&\; f \qquad \text{in the domain $\Omega$}\\
    \left(\vec{\nabla} \Phi\right) \cdot \hat{n} \;=&\; g \qquad \text{on the domain boundary $\partial\Omega$}
\end{align}


*

*We start by integrating the first equation over the entire domain:


\begin{equation}
    \int_{\Omega} \; \nabla^2 \Phi \; dV \;=\; \int_{\Omega} \; f \; dV 
\end{equation}


*We can then apply the Divergence theorem:


\begin{align}
    \int_{\Omega} \; \nabla^2 \Phi \; dV \;=&\; \int_{\Omega} \vec{\nabla}\cdot \left(\vec{\nabla} \Phi\right) dV \\
    =&\;\int_{\partial\Omega} \left(\vec{\nabla} \Phi\right) \cdot \hat{n} \; dA
\end{align}
where $\hat{n}$ is the unit vector that is perpendicular outward from the domain. 


*Note that our boundary equation tells us $\left(\vec{\nabla} \Phi\right) \cdot \hat{n} = g$, and substitute it:


\begin{align}
    \int_{\Omega} \; \nabla^2 \Phi \; dV \;=&\; \int_{\partial\Omega} \left(\vec{\nabla} \Phi\right) \cdot \hat{n} \; dA \\
    =&\;\int_{\partial\Omega} g \;dA
\end{align}


*The above relations lead to this, the compatibility constraint:


\begin{equation}
    \int_{\partial\Omega} g \;dA \;=\; \int_{\Omega} \; f \; dV 
\end{equation}
The physical interpretation of this is that the sum of the sources/sinks of flux in the interior balances the net flux through the boundary. 
How the Compatibility Constraint applies to this problem:
\begin{align}
    f \;=\; f(r,\theta) \;=&\; \frac{A}{8 \pi r^4} \csc\theta \left(1+3 \cos(2 \theta)\right) \\ \\
    g \;=\; g(r,\theta) \;=&\; 
\begin{cases}
    g_A(r) = \vec{\nabla} \Phi \cdot \left(-\hat{\theta}\right) = -\frac{A \cos (\theta_{min})}{2 \pi r^3},  & \theta = \theta_{min} \\[2ex]
    g_B(\theta) = \vec{\nabla} \Phi \cdot \left(+\hat{r}\right) = -\frac{A \sin (\theta)}{4 \pi R_{max}^3},  & r = R_{max} \\[2ex]
    g_A(r) = \vec{\nabla} \Phi \cdot \left(+\hat{\theta}\right) = \frac{A \cos (\theta_{max})}{2 \pi r^3},  & \theta = \theta_{max} \\[2ex]
    g_D(\theta) = \vec{\nabla} \Phi \cdot \left(-\hat{r}\right) = \frac{A \sin (\theta)}{4 \pi R_{min}^3},  & r = R_{min} \\[2ex]
\end{cases}
\end{align}
Two dimensions
The problem is axisymmetric, so we could perform the integrals in 3D, but for the sake of simplicity, we'll do them in 2D. (I did both and got the same result.) 
The 2D version of the above constraint is:
\begin{equation}
    \int_{\partial\Omega} g \;ds \;=\; \int_{\Omega} \; f \; dA 
\end{equation}
We can compute the 2D domain integral of $f(r,\theta)$:
\begin{align}
    \int_{\Omega} f \;dA \;=&\; \int_{R_{min}}^{R_{max}} \int _{\theta_{min}}^{\theta_{max}} \; f(r,\theta)\; r \; dr \; d\theta\\
    =&\; \dots \\
    =&\; \frac{A \left(R_{max}^2 - R_{min}^2\right)}{8 \pi R_{max}^2 R_{min}^2} \left[ 3 \cos \theta_{max} - 3 \cos \theta_{min} + 2 \ln\left(\frac{\tan\left(\theta_{max}/2\right)}{\tan\left(\theta_{min}/2\right)}\right) \right]
\end{align} 
We can also compute the 1D boundary integral of $g(r,\theta)$ around the rectangular boundary:
\begin{align}
    \int_{\partial\Omega} g \;dA \;=&\; \int_{R_{min}}^{R_{max}} g_A(r) \; dr + \int_{\theta_{min}}^{\theta_{max}} g_B(\theta) R_{max} \; d\theta + \int_{R_{max}}^{R_{min}} g_C(r) \; dr + \int_{\theta_{max}}^{\theta_{min}} g_D(\theta) R_{min} \; d\theta\\
    =&\; \dots \\
    =&\; 0
\end{align} 
Equating the two gives:
\begin{equation}
3 \cos \theta_{max} - 3 \cos \theta_{min} + 2 \ln\left(\frac{\tan\left(\theta_{max}/2\right)}{\tan\left(\theta_{min}/2\right)}\right) = 0
\end{equation}
From this it follows that
\begin{equation}
    \cos \theta_{min} - \cos \theta_{max} = \frac{2}{3} \ln\left(\frac{\tan\left(\theta_{max}/2\right)}{\tan\left(\theta_{min}/2\right)}\right)
\end{equation}
This has no real solution, implying that the problem is not well-posed. This is unsatisfactory. 
Three dimensions
We can try this in 3D, since the problem is axisymmetric. This shouldn't make a difference, but it's worth a shot:
We can compute the 3D domain integral of $f(r,\theta)$...
\begin{align}
    \int_{\Omega} f \;dV \;=&\;  \int_0^{2 \pi}\int_{\theta_{min}}^{\theta_{max}} \int_{R_{min}}^{R_{max}} \; f(r,\theta)\; r^2 \sin\theta \; dr \; d\theta \; d\phi\\
    =&\; \dots \\
    =&\; \frac{A}{4} \left( \frac{R_{max}-R_{min}}{R_{max}R_{min}}\right)  \big( \theta_{max}-\theta_{min} + 3 \cos\left(\theta_{max}\right) \sin\left(\theta_{max}\right) - 3 \cos\left(\theta_{min}\right) \sin\left(\theta_{min}\right)\big)
\end{align} 
...And the 2D surface integral of $g(r,\theta)$:
\begin{align}
    \int_{\partial\Omega} g \;dA \;=&\; \int_0^{2 \pi} \int_{R_{min}}^{R_{max}} r \sin\theta \; g_A(r) \; dr + \int_0^{2 \pi} \int_{\theta_{min}}^{\theta_{max}} r^2 \sin\theta \; g_B(\theta) R_{max} \; d\theta + \int_0^{2 \pi} \int_{R_{max}}^{R_{min}} r \sin\theta \; g_C(r) \; dr + \int_0^{2 \pi} \int_{\theta_{max}}^{\theta_{min}} r^2 \sin\theta \; g_D(\theta) R_{min} \; d\theta\\
    =&\; \dots \\
    =&\; \frac{A}{4} \left( \frac{R_{max}-R_{min}}{R_{max}R_{min}}\right)  \big( \theta_{max}-\theta_{min} + 3 \cos\left(\theta_{max}\right) \sin\left(\theta_{max}\right) - 3 \cos\left(\theta_{min}\right) \sin\left(\theta_{min}\right)\big)
\end{align} 
By inspection, we can see that the compatibility constraint is satisfied. Thus, the 3D analysis of the compatibility constraint reveals that the problem is well-posed and has a unique solution (up to an addative constant) for any $0 \le \theta_{min} < \theta_{max} \le \pi$. 
Summary
From the 2D analysis, the compatibility constraint implies that there is no physical solution, so the problem is ill-posed.  
From the 3D analysis, (which is identical to the 2D case except for the integral over $phi$), there are infinitely many solutions ($0 \le \theta_{min} < \theta_{max} \le \pi$) which satisfy the constraint. This implies that the problem is well-posed for any physical domain. 
Physically, this should be a well-posed problem. (Since posting this question, I have other reasons to know this to be true.) I must have made a mistake in my 2D analysis with the integration weights or something. The reason the 2D case fails is because the area element for a plane of constant azimuth is $dA = r \; dr \; d\theta$, while the volume element for the 3D case is $dV = r^2 \sin\theta \; dr \; d\theta \; d\phi$. The $\sin\theta$ makes a big difference for the $\csc\theta$ in $f(r,\theta)$. 
Still, rep to anyone who can find the mistake in my analysis of the 2D case.
