# Solve Laplace equation in Cylindrical - Polar Coordinates

Hey mathematica stackexchange!! I've got a (possibly stupid) problem. I've tried many things to no avail, and I've read every post I've found on Laplace's equation.

Background: I'm trying to find the capacitance per unit length of the following system (where the image is a perpendicular view of an infinite cylinder):

My approach is to first find the potential using

$$\nabla^2_{polar}V=\frac{V^{(0,2)}(\rho ,\phi )}{\rho^2}+\frac{V^{(1,0)}(\rho ,\phi)}{\rho }+V^{(2,0)}(\rho ,\phi )=0$$

Then, I would take the (negative) gradient and find $\vec{E}$ and, using that $\sigma = \varepsilon_{0}*(\vec{E} \cdot \hat{n})$, I'd get $\sigma$ in the upper left surface and integrate it to get $Q$ (per unit length). With that, $C=\frac{Q}{\Delta V}$.

The answer is supposed to be $\frac{\varepsilon_0}{\pi}*Log[2]\sim1.9$ $pF/m$ (given in the article from where my professor adapted the problem).

Mathematica implementation:

R = 1; V0 = 1; V1 = 0; e0 = 8.854187817*^-12;
regionCyl =
ImplicitRegion[
0 <= r <= R && 0 <= p <= 2 Pi, {r, p}];

laplacianCil = Laplacian[V[r, p], {r, p}, "Polar"];
boundaryConditionCil = {DirichletCondition[
V[r, p] == V0, r == R && 0 <= p <= Pi/2],
DirichletCondition[
V[r, p] ==
V1, r == R && Pi <= p <= 3/2 Pi]};

solCyl = NDSolveValue[{laplacianCil == 0, boundaryConditionCil},
V, {r, 1*^-12, R}, {p, 0, 2 Pi}, MaxSteps -> Infinity];

electricFieldCyl = -Grad[solCyl[r, p], {r, p},
"Polar"];
sigmaCyl = (Dot[
electricFieldCyl, -{1, 0}] /. {r ->
R})*e0;
Q0Cyl = NIntegrate[sigmaCyl, {p, 0, Pi/2}];
capacitanceCyl = Abs[Q0Cyl]/Abs[V0 - V1]


With this, I get $C = 5.644674742129655*10^{-12}$.

The questions:

1. How should I go about solving this problem? Is my implementation correct? I mean, are boundary conditions and the equation correctly passed to NDSolveValue? (ignoring physics).

2. Calculation time is extremely fast (less than 2 secconds to run all the code + some plots). So much so that I think the NDSolveValue is not correctly solving my problem. How can I increase resolution?

I tried to add another Boundary Condition for $V[\infty,\phi]=0, 0 < \phi <2 \pi$ with

DirichletCondition[V[r, p] == 0, r == 1*^6*R && 0 <= p <= 2*Pi];


However, NDSolve fails to solve the equation.

(I believe I'm missing $\sigma$ on the outer side of the cylinder. I'm not sure how to evaluate this. This part of the problem is not for this forum so I tried to limit my question to Mathematica. But, by all means, correct my physics if they're wrong)

• One problem I see is that your resultant potential is not periodic in p. If you plot solCyl, you will see that the values at p = 2 Pi are not the same as at p = 0, but then your bc's are not requiring that. MMa does seem to be honoring the bc's you have specified, however. Your V at your minimum r should be pretty constant with p, but is not, but haven't specified that bc either. As an aside, the capacitance should only depend on the geometry of the system and not the applied potentials, so if C is all you want, you could try more symmetric values as a check. – Bill Watts May 31 '18 at 21:19
• @BillWatts "The values at p=2 Pi are not the same as at p=0" : see my Warning at the end of this answer – andre314 Jun 1 '18 at 0:51
• I'm pretty sure the answer has to be single valued in polar angle from a physics point of view. MMa doesn't care though if not specified. – Bill Watts Jun 1 '18 at 1:10
• @andre: I've edited my post with the exact (analytic) result that I was given. I'll be waiting for more info on your progress!! Thanks! – Peanut14 Jun 1 '18 at 1:50
• @BillWatts I tried adding a DirichletCondition for V[R, 2Pi] ==V0, but it didn't change my answer. How would you solve the problem if you had total freedom to do so? BTW, I'm trying to solve the problem with V1=-1; V0=1 for that same reason. But Mma converges on the same answer, no matter what potentials I select (which is what should happen, I believe). – Peanut14 Jun 1 '18 at 1:56

I think it is better to use Cartesian coordinates because one does not have to deal with the periodicity in p. To have control on the meshing of the the region, we tell Mathematica explicitly to discretize it. PrecisionGoal -> 6 controls the meshing at the boundary. But this does not always work. Alternatively, one can use MaxCellMeasure -> 0.001 or a MeshRefinementFunction as done in the calculation below. The MeshRegion (FullForm[regionCyl]) is then passed directly to NDSolveValue. user21 pointed out that one may get a higher quality solution when using ElementMesh because it uses second-order mesh elements (reference). To specify the angle in the boundary conditions using {x,y}-coordinates you need ArcTan with two arguments (reference). ArcTan[y/x] would only cover an interval from -Pi/2 to Pi/2. The electrostatic potential (sol) as a function of x and y is calculated by NDSolveValue. The electric field is the negative gradient of the electrostatic potential. The plots below visualize the electric field lines together with the potential. The charge in a region of interest is given by the scalar product of its normal vectors and the electric field (sigmaCyl) integrated along the (closed) boundary of the region. This flux normal to the boundary (sigmaCyl) is plotted below from -2 Pi to 2 Pi. Integration of the flux is carried out for one of the plates (from 0 to p0). The electrical field exactly on the boundary is not completely covered by the mesh due to numerical inaccuracies. That is why I use the field close to the boundary at 0.999 R. The capacitance (capacitanceCyl) of the structure is given by its charge per voltage.

Clear[sigmaCyl]
R = 1; V0 = 1; V1 = 0; e0 = 8.854187817*^-12; p0 = Pi/2;
regionCyl = DiscretizeRegion[ImplicitRegion[Sqrt[x^2 + y^2] <= R, {x, y}], PrecisionGoal -> 6]
laplacian = Laplacian[V[x, y], {x, y}];
boundaryCondition = {
DirichletCondition[V[x, y] == V0, 0 < ArcTan[x, y] < p0],
DirichletCondition[V[x, y] == V1, -Pi < ArcTan[x, y] < -Pi + p0]};
sol = NDSolveValue[{laplacian == 0, boundaryCondition}, V, {x, y} \[Element] regionCyl];
electricField[x_, y_] = -Grad[sol[x, y], {x, y}];
Row[{Show[
DensityPlot[sol[x, y], {x, y} \[Element] regionCyl, ColorFunction -> "TemperatureMap", ImageSize -> Medium],
StreamPlot[electricField[x, y], {x, y} \[Element] regionCyl, StreamStyle -> Black]],
Plot3D[sol[x, y], {x, y} \[Element] regionCyl, ColorFunction -> "TemperatureMap", BoxRatios -> {1,1,1}, ImageSize -> Medium]}]
sigmaCyl[p_] = electricField[0.999 R Cos[p], 0.999 R Sin[p]].{Cos[p], Sin[p]}*e0;
Plot[sigmaCyl[p], {p, -2 Pi, 2 Pi}]
Q0Cyl = NIntegrate[sigmaCyl[p], {p, 0, p0}, AccuracyGoal -> 5];
capacitanceCyl = Abs[Q0Cyl]/Abs[V0 - V1]


Edit: As pointed out by Peanut14, to get a physically meaningful capacitance one has to consider the electric field outside of the cylinder, too. Here, a MeshRefinementFunction is used to get smaller mesh elements for r < 3. The function gets two parameters from DiscretizeRegion. The first one is a list of the coordinates of the 3 edges of each element. The second one is its area. DiscretizeRegion expects a Boolean result telling it whether it should refine the element or not. For reasons of speed this function is compiled. You may also pass a not compiled function. Then DiscretizeRegion will compile it for you. But the problem is that it does not throw an error message if it fails (as of version 11.3). Instead it just ignores the refinement function.

Clear[sigmaCyl]
R = 1; V0 = 1; V1 = 0; e0 = 8.854187817*^-12; p0 = Pi/2;
ra = DiscretizeRegion[ImplicitRegion[Sqrt[x^2 + y^2] <= 10, {x, y}]];
rb = RegionUnion[
DiscretizeRegion[ParametricRegion[r {Cos[p], Sin[p]}, {{r, 1, 1.1}, {p, 0, p0}}]],
DiscretizeRegion[ParametricRegion[r {Cos[p], Sin[p]}, {{r, 1, 1.1}, {p, -Pi, -Pi + p0}}]]];
mrf = Compile[{{vertices, _Real, 2}, {area, _Real, 0}}, If[area > 10^-2 && Norm[Mean[vertices]] < 3, True, False]];
regionCyl = DiscretizeRegion[RegionDifference[ra, rb], MeshRefinementFunction -> mrf]
laplacian = Laplacian[V[x, y], {x, y}];
boundaryCondition = {
DirichletCondition[V[x, y] == V0, 0 <= ArcTan[x, y] <= p0 && Norm[{x, y}] < 1.5],
DirichletCondition[V[x, y] == V1, -Pi <= ArcTan[x, y] <= -Pi + p0 && Norm[{x, y}] < 1.5]};
sol = NDSolveValue[{laplacian == 0, boundaryCondition}, V, {x, y} \[Element] regionCyl];
electricField[x_, y_] = -Grad[sol[x, y], {x, y}];
s[t_] = {1/Sqrt[2], 1/Sqrt[2]} + RotationMatrix[Pi/4].{Cos[t], 1.5 Sin[t]};
n[t_] = FrenetSerretSystem[s[t], t][[2, 2]](*normals to s[t]*);
Row[{Show[
DensityPlot[sol[x, y], {x, -2, 2}, {y, -2, 2}, ColorFunction -> "TemperatureMap", ImageSize -> Medium],
StreamPlot[electricField[x, y], {x, -2, 2}, {y, -2, 2}, StreamStyle -> Black],
ParametricPlot[s[t], {t, 0, 2 Pi}]],
Plot3D[sol[x, y], {x, y} \[Element] RegionIntersection[regionCyl, DiscretizeRegion[Rectangle[{-2, -2}, {2, 2}]]], ColorFunction -> "TemperatureMap", BoxRatios -> {1,1,1}, ImageSize -> Medium]}]
sigmaCyl[t_] = n[t].electricField @@ s[t]*e0;
Plot[sigmaCyl[t], {t, 0, 2 Pi}]
Q0Cyl = NIntegrate[sigmaCyl[t], {t, 0, 2 Pi}, AccuracyGoal -> 5];
capacitanceCyl = Abs[Q0Cyl]/Abs[V0 - V1]


Now, that is the capacitance.

Edit: The situation your professor is referring to is slightly different. A closed cylinder is divided in 4 segments by two perpendicular planes. One is along the symmetry axis the other with variable position (adjustable by p0). The upper-right segment is on the potential V0. The other 3 segments are on ground. Now, the charge on the lower-left segment is always the same for a given voltage independent of the position of the second plane. The so-called cross-capacitance between the lower-left and the upper-right segment is ε0/π∗Log[2]. This holds even for arbitrarily shaped cross-sections as long as they are mirror symmetric. Such a configuration is believed to result in a very stable capacitor. The original paper of Thompson and Lampard is not freely accessible but here is an open access paper that explains a bit of the context.

Clear[sigmaCyl]
R = 1; V0 = 1; V1 = 0; e0 = 8.854187817*^-12; p0 = 0.5 Pi /2;
regionCyl = DiscretizeRegion[ImplicitRegion[Sqrt[x^2 + y^2] <= R, {x, y}], PrecisionGoal -> 6];
laplacian = Laplacian[V[x, y], {x, y}];
boundaryCondition = {DirichletCondition[V[x, y] == V0, 0 < ArcTan[x, y] < p0], DirichletCondition[V[x, y] == V1, True]};
sol = NDSolveValue[{laplacian == 0, boundaryCondition}, V, {x, y} \[Element] regionCyl];
electricField[x_, y_] = -Grad[sol[x, y], {x, y}];
Row[{Show[
DensityPlot[sol[x, y], {x, y} \[Element] regionCyl, ColorFunction -> "TemperatureMap", ImageSize -> Medium, PlotRange -> All],
StreamPlot[electricField[x, y], {x, -1, 1}, {y, -1, 1}, StreamStyle -> Black],
Graphics[{Thick, Circle[{0, 0}, 1, {0, p0}], Circle[{0, 0}, 1, {-Pi, -p0}], Dashed, Line[{{Cos[p0], -1}, {Cos[p0], 1}}], Line[{{-1, 0}, {1, 0}}]}]],
Plot3D[sol[x, y], {x, y} \[Element] regionCyl, ColorFunction -> "TemperatureMap", BoxRatios -> {1, 1, 1}, ImageSize -> Medium, PlotRange -> All]}]
sigmaCyl[p_] = electricField[0.9999 R Cos[p], 0.9999 R Sin[p]].{Cos[p], Sin[p]}*e0;
Plot[sigmaCyl[p], {p, -2 Pi, 2 Pi}]
Q0Cyl = NIntegrate[sigmaCyl[p], {p, -Pi, -p0}, AccuracyGoal -> 5];
capacitanceCyl = Abs[Q0Cyl]/Abs[V0 - V1]


Note that the vertical line does not have to be in the middle. The cross-capacitance is always the same. In general, the electric field between the lower-left and the upper-right segment outside of the cylinder has to be considered, too. But it is smaller due to the other segments. In practical situations the segmented cylinder is shielded with a non-segmented cylinder around it that is on ground potential.

• @andre ah, yes. The value I got is close to the one that you calculated. That was for a plate capacitor with plates of length R Sqrt[2] and distance R Sqrt[2], right? That seems to make sense. – Matthias Bernien Jun 2 '18 at 12:14
• Thanks for the answer! Yesterday I went to his office, but he was not there. I'll have to wait until Monday to talk this through with him. I too suspect of the "correct" answer, but was unsure about my &5.6 pF/m& answer as well. Question: by evaluating sigmaCyl at 0.999*R, are you considering a plaque of thickness 0.001? If so, wouldn't I need to also solve for the $\sigma$ on the outer face of the plaque? – Peanut14 Jun 2 '18 at 13:31
• BTW, I'm using MMa 10 and I get this error Range specification {x,y}\[Element]regionCyl is not of the form {x, xmin, xmax} (from Stream plot). I don't really need a solution for that (I can work around it). Just thought I'd let you know :) . – Peanut14 Jun 2 '18 at 13:35
• @Peanut14 You are welcome! The electrical field exactly on the boundary is not completely covered by the mesh due to numerical inaccuracies. That is why I used 0.999 R. The outer face is actually not existent in the model because the potential is only calculated within the region defined by regionCyl. – Matthias Bernien Jun 2 '18 at 15:11
• @MatthiasBernien I understand the model doesn't cover that face. But, should it? I mean, this is now phyisics. But, how would I extend $V$ to cover all of space (or at least the interior of the cylinder and a region around it). Which version of MMa are you using? – Peanut14 Jun 2 '18 at 16:39

Here is a analysis of all the problems related to Mathematica in your question.

Briefly summarized, These are 3 problems :

• the lack of build-in tools to visualize the potential and vector field in polar coordinates.

• boundary problems : Whatever is the real geometry you are interested (not clear in your question, especially you want to optain 1.9 pF/m), there are boundaries that are not expected (compared to your description of the geometry). This will become clear once we will have the tools to visualize the vector field

• There is also a difficulty due with the fact that Grad[potential] return a pair of Interpolation functions and not a unique interpolation function that returns a pair of values.

Visualisation tools

You code (exactly) :

R = 1; V0 = 1; V1 = 0; e0 = 8.854187817*^-12;
regionCyl =
ImplicitRegion[
0 <= r <= R && 0 <= p <= 2 Pi, {r, p}];

laplacianCil = Laplacian[V[r, p], {r, p}, "Polar"];
boundaryConditionCil = {DirichletCondition[
V[r, p] == V0, r == R && 0 <= p <= Pi/2],
DirichletCondition[
V[r, p] ==
V1, r == R && Pi <= p <= 3/2 Pi]};

solCyl = NDSolveValue[{laplacianCil == 0, boundaryConditionCil},
V, {r, 1*^-12, R}, {p, 0, 2 Pi}, MaxSteps -> Infinity];

electricFieldCyl = -Grad[solCyl[r, p], {r, p},
"Polar"];
sigmaCyl = (Dot[
electricFieldCyl, -{1, 0}] /. {r ->
R})*e0;
Q0Cyl = NIntegrate[sigmaCyl, {p, 0, Pi/2}];
capacitanceCyl = Abs[Q0Cyl]/Abs[V0 - V1]


The potential :

potentialSquareRepresentation=ContourPlot[solCyl[r, p], {r,p} \[Element] solCyl["ElementMesh"]
, ColorFunction -> "Temperature"
,Contours-> 20
, PlotLegends -> Automatic
];
potentialCylindricalRepresentation=Show[
potentialSquareRepresentation /. GraphicsComplex[array1_, rest___] :>
GraphicsComplex[(#[[1]] {Cos[#[[2]]],Sin[#[[2]]]})& /@ array1, rest],
PlotRange -> Automatic
]


The field : , thanks to Matthias

electricField1[r_, p_] = -Grad[solCyl[r, p ], {r, p}, "Polar"];
electricField2[x_, y_] = TransformedField["Polar" -> "Cartesian", electricField1[r, p + Pi], {r, p } -> {x, y}] /. ArcTan[x_,y_]:> ArcTan[-x,-y];
fieldCylindricalRepresentation=StreamPlot[electricField2[x, y], {x, -1, 1}, {y, -1, 1}, StreamStyle -> Black]


The mesh (the one NDSolve really used) :

meshSquareRepresentation= solCyl["ElementMesh"]["Wireframe"];
meshCylindricalRepresentation=Show[meshSquareRepresentation /. GraphicsComplex[array1_, rest___] :>
GraphicsComplex[(#[[1]] {Cos[#[[2]]],Sin[#[[2]]]})& /@ array1, rest],
PlotRange ->  {{-1,1},{-1,1}}
]


a example of superposition :

Show[potentialCylindricalRepresentation,fieldCylindricalRepresentation]


Boundary problems

As one can see on the graphics, the boundaries are :

• The two expected quarters of circle , the lower at 0 Volts, the upper at 1 Volt. That's fine.
• The boundary p=0 (ie angle=0). This boundary is a problem. It is not specified. In that case NDSolve takes the default boundary condition : Neuman=0, that is to say no field transverse to the boundary. This is clearly visible when one observes the field lines.
• There are also the two quarters of circle complementary to the two quarters that have been specified. Mathematica has seen a boundary because it is the limit of the domain. Once again the default boundary condition Neumann=0 has been used (see the field lines)

... To be continued ...

• This is certainly not the mesh NDSolve used. Have a look at: solCyl["ElementMesh"]["Wireframe"] – user21 Jun 4 '18 at 6:37
• @user21 I have corrected the mesh. Thanks. – andre314 Jun 4 '18 at 6:52
• Working still in progress. At the moment, I have put aside the eventual problem at the center (by putting a real hole of diameter 0.4 ) and I'm trying to make the boundary p=0 disappear, ie impose continuity on the potential and the field at this boundary. Not a simple affair, unless I miss something trivial. – andre314 Jun 5 '18 at 15:56
• Until now, I didn't succeed in applying Periodic Neuman condition (for the continuity of the field). I know it works in 1D, but In 2D ???. Furthermore the aim is to impose also continuity of potential, that is to say to use PeriodicBoundaryCondition[] – andre314 Jun 5 '18 at 16:01
• It turns out that the use of PeriodicBoundaryCondition[] gives the choice to choose the direction of propagation of the solution (see the doc of PeriodicBoundaryCondition[] item "Possible issue"). One can calculate the two corresponding solutions. I wonder If the solution that respect continuity of the Field (normaly imposed by Neuman) is no more than a linear combinaison of these two solutions. Suspens... – andre314 Jun 5 '18 at 16:07

The last geometry we can solve the problem analytically. We set V = V0 for p from 0 to Pi/2, and V = V1 from Pi/2 to 2Pi.

Clear["Global*"]

pde = Laplacian[V[r, p], {r, p}, "Polar"] == 0;


Separate Variables

V[r_, p_] = R[r] P[p];

Expand[(r^2*pde)/V[r, p]]

P''[p]/P[p] + (r^2 R''[r])/R[r] + (r R'[r])/R[r] == 0


Each section must be equal to a constant. We know the solution must be periodic in p so choose

peq = P''[p]/P[p] == -a^2;

DSolve[peq, P[p], p] // Flatten
{P[p] -> C[2]*Sin[a*p] + C[1]*Cos[a*p]}

p1 = P[p] /. % /. {C[1] -> c1, C[2] -> c2}


The r equation becomes

req = -a^2 + (r^2 R''[r])/R[r] + (r R'[r])/R[r] == 0;

DSolve[req, R[r], r] // Flatten // TrigToExp;

r1 = R[r] /. % /. {C[1] -> c3, C[2] -> c4}

r1 // Collect[#, r^_] &
(*(c3/2 - (I*c4)/2)/r^a + r^a*(c3/2 + (I*c4)/2)*)

r1 = % /. {c3/2 - (I*c4)/2 -> c3, c3/2 + (I*c4)/2 -> c4}
(*c3 r^-a+c4 r^a*)

Vin[r_, p_] = r1 p1
(*(c3/r^a + c4*r^a)*(c1*Cos[a*p] + c2*Sin[a*p])*)


Vin is bounded at r = 0, and single valued in p, which requires

c3 = 0
c4 = 1
a = n
$Assumptions = n \[Element] Integers  We set c4 to 1 to combine it with c1 and c2. Vout is bounded at r = Infinity requiring c8 = 0 c7 = 1  We end up with Vin[r, p] (*r^n (c1 Cos[n p] + c2 Sin[n p])*) Vout[r, p] (*r^-n (c5 Cos[n p] + c6 Sin[n p])*)  Work with the solution Vin for r < R At r = R, V= V1 0 <=p <= Pi/2, and V0 otherwise Use orthogonality to match the boundary at r = R and solve for the c constants. the n=0 term eq0 = Integrate[V0, {p, 0, Pi/2}] + Integrate[V1, {p, Pi/2, 2*Pi}] == R^0*Integrate[c0, {p, 0, 2*Pi}]//FullSimplify Solve[%, c0]; c0 = c0 /. %[[1]];  eq1 mult by sin and integrate eq1 = Integrate[V0*Sin[n*p], {p, 0, Pi/2}] + Integrate[V1*Sin[n*p], {p, Pi/2, 2*Pi}] == R^n*Integrate[(c1*Cos[n*p] + c2*Sin[n*p])*Sin[n*p], {p, 0, 2*Pi}]//FullSimplify;  eq2 mult by cos and integrate eq2 = Integrate[V0*Cos[n*p], {p, 0, Pi/2}] + Integrate[V1*Cos[n*p], {p, Pi/2, 2*Pi}] == R^n*Integrate[(c1*Cos[n*p] + c2*Sin[n*p])*Cos[n*p], {p, 0, 2*Pi}]//FullSimplify; Solve[eq1, c2] // Flatten // FullSimplify; c2 = c2 /. %; Solve[eq2, c1] // Flatten // FullSimplify; c1 = c1 /. %;  Put in some values R = 1 V0 = 1 V1 = 0 Vin[r, p] // FullSimplify (*(2 r^n Sin[(Pi n)/4] Cos[n (p - Pi/4)])/(Pi n)*)  The full solution is the c0 term plus the sum of the above over integer n. c0 (*1/4*)$Assumptions = r >= 0 && p \[Element] Reals

Vin[r_, p_] = 1/4 + (2/Pi)*Sum[(r^n*Sin[(Pi*n)/4]*Cos[n*(p - Pi/4)])/n, {n, 1, Infinity}]//FullSimplify
(*-((I*(2*Log[1 - r/E^(I*p)] - 2*Log[1 - (I*r)/E^(I*p)] - 2*Log[1 - E^(I*p)*r] + 2*Log[1 + I*E^(I*p)*r] + I*Pi))/
(4*Pi))*)


MMa successfully finds a closed form solution to the infinite sum. It looks very complex, but plotting shows that it is a real expression.

Electric field in r direction

Efrin[r_, p_] = -D[Vin[r, p], r] // FullSimplify
(*-((I*E^(2*I*p)*((r^2 + 1)*Sin[p] + (r^2 + 1)*Cos[p] - 2*r))/(Pi*(-r + E^(I*p))*(E^(I*p) - I*r)*(-1 + E^(I*p)*r)*
(E^(I*p)*r - I)))*)

Plot[Efrin[R, p], {p, 0, 2 Pi}]


Charge density

Sigma[p_] = -e0 (Efrin[R, p] // FullSimplify)
(*-(e0/(-(Pi*Sin[p]) - Pi*Cos[p] + Pi))*)


Surprisingly simple expression for the charge density. Calculate the total q/length for the section opposite the potential V1.

q = Integrate[Sigma[p], {p, Pi, (3*Pi)/2}]

-((e0*Log[2])/Pi)

Cap = Abs[q/(V1 - V0)]
(*(e0*Log[2])/Pi*)
`

What is interesting is that the integral for Sigma over the p limits of V1 does not converge.