# Laplace PDE in a polar coordinate system

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

• A previous version of this question was closed. To avoid the same fate with this version of the question, you should show what your specific issue in Mathematica is and how you have so far tried to address it.
– Jens
Commented Mar 30, 2015 at 16:47
• @Jens I corrected the question Commented Mar 30, 2015 at 18:08
• I think you are looking for the last equation on the slide you shared. Commented Mar 30, 2015 at 18:32
• This PDF shows what to do at $r=0$ to avoid blow up, which is to set one of the constant of integrations to zero. May be you can follow this derivation and apply your numerical solution to the resulting solution. Commented Mar 31, 2015 at 14:20

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];