Laplace equation in a rectangle, Dirichlet to Neumann map

Question

Consider the problem

$$\Delta u = 0,(x,y) \in \Omega=(0,1)^2 \\ u(0,y)=u(1,y)=1 \\ u(x,0)=u(x,1)=0.$$

This problem can be solved exactly. The solution is

$$u(x,y)=4\sum_{\text{odd } n} \frac{1}{n\pi} \sin(n\pi x) \left ( 1 - \frac{\sinh(n\pi y)+\sinh(n\pi(1-y))}{\sinh(n\pi)} \right ).$$

Note that the boundary conditions on the left and right edges are not strictly speaking satisfied; these are only satisfied in a weak sense. Essentially the right way to think about this is to round the corners out into quarter-circles of a small radius $r$, make the boundary value go smoothly (but rapidly) from $0$ to $1$ around each corner, and then take the limit $r \to 0$. In the limit $r \to 0$ there is a Gibbs phenomenon appearing which spoils the boundary conditions.

I am trying to analytically compute or approximate $\frac{\partial u}{\partial y}(x,1)$. Unfortunately, it appears that the convergence properties of this sum are too bad to allow term-by-term differentiation with respect to $y$ (which is not a huge surprise, because the boundary conditions are singular at the corners). How can this be circumvented to extract the boundary derivative data?

Answer 1

If we naively compute the derivative term by term we end up with a divergent series. However the partial sums of this series will oscillate about the true solution. We can obtain an approximation by adding an exponential regulator which leads to the following approximation

$$\frac{\partial u}{\partial y}(x,1) \approx -4\sum_{n\text{ odd}}\sin(n\pi x) \tanh(\pi n/ 2)e^{-\epsilon n}$$

To check how good this approximation is we need the numerical solution. Luckily this is very easy to obtain in Mathematica:

(* Solve the PDE numerically *)
equation = Laplacian[u[x, y], {x, y}] == 0;
dirichletconditons = {
   DirichletCondition[u[x, y] == 0, y == 0], DirichletCondition[u[x, y] == 0, y == 1],
   DirichletCondition[u[x, y] == 1, x == 0], DirichletCondition[u[x, y] == 1, x == 1]};
domain = Rectangle[{0, 0}, {1, 1}];
method = {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}
f = NDSolveValue[{equation, dirichletconditons}, u, {x, y} \[Element] domain, Method -> method];
dfdy = Derivative[0, 1][f];
Plot3D[f[x, y], {x, 0, 1}, {y, 0, 1}]

For the approximation we fix $\epsilon$ and sum up to high enough $n_{\rm max}$ such that higher order terms are negligible

(* Compare numerical solution to analytical approximation *)
\[Epsilon] = 0.01;
nmax = Floor[5.0/\[Epsilon]];
dfdyapprox[x_] := Sum[-4 Sin[\[Pi] n x] Tanh[\[Pi] n/2] Exp[-\[Epsilon] n], {k, 1, nmax, 2}];
Plot[{dfdy[x, 1], dfdyapprox[x]}, {x, 0, 1}]

which gives an almost perfect match (the two curves lie on top of each other here)

$~~~~~~~~~~~~~$

Now if we make the approximation $\tanh(x)\approx 1$ in the sum above we can evaluate the sum and take $\epsilon \to 0$ to get $$\frac{\partial u}{\partial y}(x,1) \approx \frac{4 \sin (\pi x)}{\cos (2 \pi x)-1}$$ which turns out to be a quite good approximation

and this can be improved even further by adding $-4\sum_{n=1,3,5,(2N+1)}\sin(n\pi x)[\tanh(n\pi/2)-1]$ for some small $N$. Turns out adding just one term gives a result that accurate to $0.1\%$ for most $[0,1]$:

$$\frac{\partial u}{\partial y}(x,1) \approx 4 \sin (\pi x) \left(\frac{1}{\cos (2 \pi x)-1}+1-\tanh \left(\frac{\pi }{2}\right)\right)$$

Answer 2

The derivative series is just fine for $0 < y < 1$: $$ u_{y}(x,y) = 4\sum_{\mbox{odd $n$}}\sin(n\pi x)\frac{\cosh(n\pi(1-y))-\cosh(n\pi y)}{\sinh(n\pi)} $$ For $\delta < y < 1-\delta$ the fraction in the sum is bounded by $C e^{-n\pi\delta}$.