# How can I solve the 2D Laplace equation with Neumann boundary conditions? [closed]

I would like to solve the following boundary value problem for $u(x,y)$:

### PDE

$\quad \quad \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0$

### Domain

$\quad \quad 0 \le x \le L$

$\quad \quad 0 \le y \le H$

### Boundary conditions

$\quad \quad \frac{\partial u}{\partial x} \big|_{x=0} =\; B$

$\quad \quad \frac{\partial u}{\partial x} \big|_{x=L} =\; B$

$\quad \quad \frac{\partial u}{\partial y} \big|_{y=0} =\; 0$

$\quad \quad \frac{\partial u}{\partial y} \big|_{y=H} =\; 0$

$B$ is a constant.

FYI: Unless I'm mistaken, the general solution for this boundary value problem (BVP) takes the form:

\begin{equation} u(x,y) = - \sum_{n=0}^\infty C_n \left( \frac{B}{\lambda_n}\right) \left( \frac{\sinh \left(\frac{\lambda_n L}{2}\right)}{\sinh \left(\lambda_n L\right)}\right) \cos\left( \lambda_n y \right) \left[ 2 \sinh\left(\frac{\lambda_n L}{4}\right) e^{\lambda_n\left(\frac{L}{4} + x\right)} + e^{\lambda_n\left(\frac{L}{2} + x\right)} \right], \end{equation}

where $\lambda_n \equiv \frac{n \pi}{H}$, and $C_n$ are undetermined real coefficients.

I can solve this by hand easily enough via separation of variables, but I suspect Mathematica ought to be able to solve this.

### Code

(* Start with a clean kernel *)
ClearAll["Global*"];
Remove["Global*"];

(* The PDE *)
LaplaceEqn = D[u[x, y], x, x] + D[u[x, y], y, y] == 0;

(*The Boundary Conditions*)
minXBC = Derivative[1, 0][u][0, y] == b;
maxXBC = Derivative[1, 0][u][L, y] == b;
minYBC = Derivative[0, 1][u][x, 0] == 0;
maxYBC = Derivative[0, 1][u][x, H] == 0;

soln = DSolve[{LaplaceEqn, minXBC, maxXBC, minYBC, maxYBC}, u, {x, y}]


Unfortunately, that code produces the following output (Mathematica 10):

 DSolve[
{(u^(0,2))[x,y]+(u^(2,0))[x,y] == 0,
(u^(1,0))[0,y] == b, (u^(1,0))[L,y] == b,
(u^(0,1))[x,0] == 0, (u^(0,1))[x,H] == 0},
u, {x, y}]


This kind of PDE is supposedly solveable by Mathematica, according to the documentation of DSolve

I can invoke DSolve in the following way...

soln = DSolve[LaplaceEqn, u, {x,y}]


...but it gives useless results:

{{u -> Function[{x, y}, C[I x  y] + C[-I x + y]]}}


I understand that that might mean something to somebody, but I want a real-valued solution in closed-form (a series expansion with undetermined weights is okay).

What am I doing wrong? How can Mathematica solve this analytically in a way that takes these boundary conditions into account?

• @march I'm guessing he copied that from the output - so that wasn't an actual error in the input... – Jens Jul 23 '15 at 17:27
• The general result for DSolve[LaplaceEqn,u,{x,y}] is indeed correct, because any arbitrary (holomorphic) functions of $x\pm i y$ satisfy the equation. But with boundary conditions, there is no closed form solution, and the series solution is not what you'll get unless you specifically ask for an expansion in terms of eigenfunctions (but that has to be put in by hand). – Jens Jul 23 '15 at 17:29
• Is your goal to obtain specifically the series form as you wrote it, and determine the coefficients? Or do you want a purely numerical solution? – Jens Jul 23 '15 at 17:33
• FYI, your original code had all of the boundary derivatives equal to zero (in which case the only solution is a constant.) I've modified the code to match the equations. – Michael Seifert Jul 23 '15 at 17:47
• @MichaelSeifert It looks like your edit of the code was the opposite of what you wanted. Anyway, to answer the question we need to know if the goal is to obtain the (exact but infinite) expansion or a numerical result. – Jens Jul 23 '15 at 17:59