Laplace PDE in a polar coordinate system

Question

I want to solve a Laplace PDE in a polar coordinate system with finite difference method, but I have a problem with boundary conditions at r = 0.

Here is what I found on the Internet: http://homepages.see.leeds.ac.uk/~amt6xw/Distance%20Learning/CFD5030/node10.html

The formula for discretization is here:

I started off from a Cartesian coordinate system, from a rectangle grid. I tried to transfrom it into polar coordinates, but I don't know how to add/define the boundary conditions at r = 0, but in r = R = 0.

R = 0.04;
n = 10;
m = 10;
Δr = R/n;
Δθ = R/m;
ϕ[n, j_] = 0;(*???*)
ϕ[i_, m] = 0;(*???*)
vars = Flatten[Table[ϕ[i, j], {i, 1, n - 1}, {j, 1, m - 1}]];
eqns = Flatten[Table[1/Δr^2 (1 - 1/(2 m)) ϕ[i - 1, j] - 
  2 ϕ[i, j] + (1 - 1/(2 m)) ϕ[i + 1, j] + 
  1/(m Δθ)^2 (ϕ[i, j + 1] - 
     2 ϕ[i, j] + ϕ[i, j - 1]) == 0, {i, 1, n - 1}, {j, 1, m - 1}]];
sol = Solve[eqns, vars][[1]];
ϕsol = 
  Interpolation[
    Flatten[Table[{i Δr, j Δθ, ϕ[i, j]}, {i, 0, n}, {j, 0, m}] /. sol, 1]];

When I solve analitycal the result is here:

But I don't know how to solve numerically

Answer 1

OK, let me kill this unanswered question. As already pointed out in the comment above, one possible treatment for $\nabla^2 U$ at $r=0$ is

$$\nabla^2 U = \frac{4 \left( u_m - u_0 \right)}{\left( \Delta r \right)^2}$$

where $u_m$ is the mean value of $u$ along $r = \Delta r$. Aside from this, there're dozens of simple mistakes in the original code sample and I'd like not to point out them one by one. The last illustration isn't correct either, it's probably a solution of Possion equation (To be more specific, the equation in another question of OP, I guess ). Anyway, the following is the fixed code. The b.c. chosen here is $ϕ(R, θ) = \sin(2 \pi θ)$.

Clear[ϕ]
R = 4/100;
n = 25;
m = 25;
Δr = R/n;
Δθ = 2 Pi/m;

ϕ[0, j_] = ϕ[0, 0];
ϕ[n, j_] = Sin[2 Pi j/m];
ϕ[i_, j_?(! 0 <= # < m &)] = ϕ[i, Mod[j, m]];

vars = Union@Flatten[Table[ϕ[i, j], {i, 0, n - 1}, {j, 0, m}]];
eqns = Flatten@
   Table[((ϕ[i - 1, j] - 2 ϕ[i, j] + ϕ[i + 1, j])/Δr^2 + 
       1/(i Δr) (ϕ[i + 1, j] - ϕ[i - 1, j])/(2 Δr) + 
       1/(i Δr)^2 (ϕ[i, j + 1] - 2 ϕ[i, j] + ϕ[i, j - 1])/ (Δθ^2)) == 0, 
     {i, 1, n - 1}, {j, 1, m}];

eqadditional = With[{ϕm = Mean@Table[ϕ[1, j], {j, 1, m}]}, (4 (ϕm - ϕ[0, 0]))/Δr^2 == 0];

solrule = NSolve[Flatten@{eqns, eqadditional}, vars][[1]];
solcoord = solrule /. (_[i_, j_] -> val_) :> {i Δr Cos[j Δθ], i Δr Sin[j Δθ], val};
bccoord = Table[{n Δr Cos[j Δθ], n Δr Sin[j Δθ], N@ϕ[n, j]}, {j, 1, m}];

ListPointPlot3D[solcoord~Join~bccoord]

Notice you can select the formula and press Ctrl+Shift+T to make it more readable.