# NDSolve returns asymmetric solution to a symmetric equation

I am trying to find the potential of a conducting cylindrical electrode by solving the Laplace's equation.

Both, the boundary conditions and the equation are symmetric w.r.t. the change $$r\to-r$$. However, the solution is not. I guess that the discrepancy is due to lack of accuracy in numerical calculations.

So the question is: how can the accuracy be improved in order to produce a symmetric outcome?

My code is given below.

Initiating the variables (decided not to change my specific numbers)

rElectrode = 0.02;
lElectrode = 0.6;
lCamera = 1.5;
rCamera = 1.5/2;
voltage = 2.*10^4;


Defining the domain for variables $$(r,z)$$

domain = RegionDifference[
Rectangle[{-rCamera, -lCamera/2}, {rCamera, lCamera/2}],
Rectangle[{-rElectrode, -lElectrode/2},
{rElectrode, lElectrode/2}]];


Constructing the solution using NDSolve

phi = NDSolveValue[{D[r*D[\[Phi][r, z], r], r] +
r*D[\[Phi][r, z], z, z] == 0,
DirichletCondition[\[Phi][r, z] == voltage,
(Abs[r] == rElectrode && -lElectrode/2 <= z <=
lElectrode/2) || (-rElectrode <= r <= rElectrode &&
Abs[z] == lElectrode/2)],
DirichletCondition[\[Phi][r, z] == 0.,
Abs[r] == rCamera || Abs[z] == lCamera/2]}, \[Phi],
Element[{r, z}, domain], "ExtrapolationHandler" -> {Indeterminate & }];


And, finally, plotting the solution

 DensityPlot[phi[r, z], Element[{r, z}, domain],
PlotRange->All,ColorFunction->"Rainbow"]


The problem somewhat persists. As mentioned in the comments, one can force more points to be used for plotting with the command PlotPoints to get smoother result. However, I still think that there are some issues related to accuracy of numerical computations.

While the solution itself looks smooth enough, this is not true for its derivatives. For instance, let us

Plot[-Derivative[1, 0][phi][r, 0], {r, -rCamera, rCamera},
PlotRange -> All,
AxesOrigin -> {0, 0}, PlotPoints -> 1000]


producing

The solution is not ideally symmetric (one can see the spike at the left leaf near $$r=-0.2$$ which is not present on the right). But this is just a minor problem in the case at hand. It looks like very few points is used by Mathematica to plot the graph near the origin. Increasing number of PlotPoints does not change the picture. So it seems to me that there is not enough accuracy present in the solution to be differentiated. I've tried to play with the AccuracyGoal,WorkingPrecision etc., but nothing worked. However, I have very little experience in numerical computing with Mathematica and may missed something big. So my question basically remains

• Are my problems accuracy-related? If so

How do I control the precision of the solution in order to

• See the symmetry to an arbitrary desirable extent
• Find smooth enough derivatives of the solution?
• If it is the plot that looks asymmetrical to you, you can improve that dramatically simply by requesting that more points be evaluated: e.g. add PlotPoints -> 100 as an option to your density plot. – MarcoB Jun 22 '15 at 19:39
• @MarcoB, that definitely solves my problem. You could maybe turn this into a short answer? – Weather Report Jun 22 '15 at 19:46
• Also, in general, if you KNOW about symmetries of the solution, you can solve for a smaller region (in this case, a half space). This is also more efficient numerically. – yohbs Jun 22 '15 at 21:50

It seems increasing the number of requested points and the accuracy of NDSolve does the trick?

rElectrode = 0.02;
lElectrode = 0.6;
lCamera = 1.5;
rCamera = 1.5/2;
voltage = 2.*10^4;
domain = RegionDifference[
Rectangle[{-rCamera, -lCamera/2}, {rCamera, lCamera/2}],
Rectangle[{-rElectrode, -lElectrode/2}, {rElectrode,
lElectrode/2}]];

phi = NDSolveValue[{D[r*D[ϕ[r, z], r], r] +
r*D[ϕ[r, z], z, z] == 0,
DirichletCondition[ϕ[r, z] ==
voltage, (Abs[r] == rElectrode && -lElectrode/2 <= z <=
lElectrode/2) || (-rElectrode <= r <= rElectrode &&
Abs[z] == lElectrode/2)],
DirichletCondition[ϕ[r, z] == 0.,
Abs[r] == rCamera || Abs[z] == lCamera/2]}, ϕ,
Element[{r, z}, domain],
"ExtrapolationHandler" -> {Indeterminate &},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 1/10000},
"IntegrationOrder" -> 8}}];

DensityPlot[phi[r, z], Element[{r, z}, domain], PlotRange -> All,
ColorFunction -> "Rainbow", PlotPoints -> 150]


ContourPlot[phi[r, z], Element[{r, z}, domain], PlotRange -> All,
ColorFunction -> "Rainbow", PlotPoints -> 150, Contours -> 30]


Plot[Table[phi[r, z], {z, -0.6, 0.6, 0.1}] // Evaluate, {r, -0.6, 0.6},
PlotPoints -> 150]


Though there are small but finite residual errors:

Plot[Table[phi[r, z] - phi[-r, z], {z, -0.6, 0.6, 0.1}] //
Evaluate, {r, -0.6, 0.6}, PlotPoints -> 150, PlotRange -> All]