How can I plot the gradient field in bipolar coordinates?

how I can plot the gradient of function V according to the attached program in bipolar coordinates. I used from ParametricPlot to show bipolar coordinates as:

Show[ParametricPlot[
Evaluate[Table[
Tooltip[{Sinh[v]/(Cosh[v] - Cos[u]), Sin[u]/(Cosh[v] - Cos[u])},
Row[{"u \[LongEqual] ", u}]], {u, 0, 2 Pi, Pi/10}]], {v, -1, 1}],
ParametricPlot[
Evaluate[Table[
Tooltip[{Sinh[v]/(Cosh[v] - Cos[u]), Sin[u]/(Cosh[v] - Cos[u])},
Row[{"v \[LongEqual] ", v}]], {v, -1, 1, 1/10}]], {u, Pi/100,
2 Pi}]]


Now I want to calculate the gradient of the function V which is defined as:

NN = 150;
V0 = 5;
tau0 = 0.5;
V = V0*Sqrt [2*(Cosh [tau] - Cos [sigma])]*(
LegendreQ[0 - 1/2, Cosh [tau0]]*LegendreP[0 - 1/2, Cosh [tau]]*
Cos [0*sigma]/LegendreP[0 - 1/2, Cosh [tau0]] +
2*Sum[LegendreQ[n - 1/2, Cosh [tau0]]*LegendreP[n - 1/2, Cosh [tau]]*
Cos [n*sigma]/LegendreP[n - 1/2, Cosh [tau0]], {n, 1, NN}])/Pi;


After taking the gradient I want to show the results in bipolar coordinate. The sample result is shown in the following.

[![enter image description here][3]][3]

• Could you check your second definition (V)? It is not written properly. Sep 10, 2020 at 11:41
• Thanks. I think it is true. Anyway, I provide again it in the original form as a picture.
– iman
Sep 10, 2020 at 12:04
• LegendreQ[n/2,x] is a complex number for x>1. I think V will also be a complex number. You can not plot complex number. Sep 10, 2020 at 13:31
• Thanks. So can we use the Legendre function of the second kind of type 3 instead? Or select domain differently to have x<=1. Please consider any correct value to plot this function
– iman
Sep 10, 2020 at 13:48
• To evaluate Q for x>1, we can use an optional argument to Mathematica's LegendreQ[] that puts the branch cut from −∞ to +1. This then makes Q real rather than complex for x>1.
– iman
Sep 10, 2020 at 18:38

Prelude: bipolar coordinates

I use $$(u, v)$$ rather than $$(\sigma, \tau)$$ throughout. The transformation is given by\begin{align*} x &= \frac{\sinh v}{\cosh v - \cos u} \\ y &= \frac{\sin u}{\cosh v - \cos u}.\end{align*}

Not so well known is the inverse transformation (which we will need for plotting):\begin{align*} u &= \tan^{-1} \frac{2y}{x^2 + y^2 - 1} \\ v &= \tanh^{-1} \frac{2x}{x^2 + y^2 + 1}.\end{align*}

The scale factors (Lamé coefficients) for both coordinates are the same:$$h_u = h_v = h = \frac{1}{\cosh v - \cos u}.$$

The local orthonormal basis is related to the standard Cartesian basis according to\begin{align*} \mathbf{a}_u &= h (-S \,\mathbf{a}_x + C \,\mathbf{a}_y) \\ \mathbf{a}_v &= h (-C \,\mathbf{a}_x - S \,\mathbf{a}_y), \\\end{align*} where \begin{align*} C &= \cos u \cosh v - 1 \\ S &= \sin u \sinh v.\end{align*}

Implementing this in Mathematica:

(* Coordinate transformations *)
xBipolar[u_, v_] := Sinh[v] / (Cosh[v] - Cos[u]);
yBipolar[u_, v_] := Sin[u] / (Cosh[v] - Cos[u]);

(* Inverse coordinate transformations *)
uBipolar[x_, y_] := ArcTan[x^2 + y^2 - 1, 2 y];
vBipolar[x_, y_] := ArcTanh[2 x / (x^2 + y^2 + 1)];
uvBipolar[x_, y_] := {uBipolar, vBipolar} @@ {x, y} // Through // Evaluate;

(* Scale factors (both are the same) *)
hBipolar[u_, v_] := 1 / (Cosh[v] - Cos[u]);

(* Abbreviations *)
cBipolar[u_, v_] := Cos[u] Cosh[v] - 1;
sBipolar[u_, v_] := Sin[u] Sinh[v];

(* Cartesian components of local orthonormal basis *)
uVectorBipolar[u_, v_] := hBipolar[u, v] {-sBipolar[u, v], cBipolar[u, v]} // Evaluate;
vVectorBipolar[u_, v_] := hBipolar[u, v] {-cBipolar[u, v], -sBipolar[u, v]} // Evaluate;


Function

I have taken $$v_0 = 1$$. We have$$\frac{V}{V_0} = \sum_{n = 0}^\infty \frac{1}{\pi} \sqrt{2 (\cosh v - \cos u)} \cdot \lambda_n \cdot \frac{Q_{n-1/2}(\cosh v_0)}{P_{n-1/2}(\cosh v_0)} \cdot P_{n-1/2}(\cosh v) \cos (n u).$$

(* lambda-bar *)
lambdaBar[0] = 1;
lambdaBar[n_] /; n > 0 = 2;

(* Expansion terms *)
v0 = 1;
term[n_][u_, v_] := (
1 / Pi
Sqrt[2 (Cosh[v] - Cos[u])]
lambdaBar[n]
LegendreQ[n - 1/2, Cosh[v0]] / LegendreP[n - 1/2, Cosh[v0]]
LegendreP[n - 1/2, Cosh[v]]
Cos[n u]
);

(* Partial sum for V/V_0 *)
partialSum[nMax_][u_, v_] := Sum[term[n][u, v], {n, 0, nMax}];


As mentioned in the comments, the LegendreQ factor isn't real, so we take the real part.

First we figure out how many terms we need to plot; I am guessing that the boundary condition which gives rise to the Fourier series is $$V/V_0 = 1$$ along the circle $$v = v_0$$:

nMaxValues = {0, 1, 2, 5};

Plot[
Table[
partialSum[nMax][u, v0] // Re
, {nMax, nMaxValues}
] // Evaluate
, {u, 0, 2 Pi}
, PlotLegends -> LineLegend[nMaxValues, LegendLabel -> "nMax"]
]


We see that 5 terms are enough. We certainly don't need 150.

(* Assume function intended for v < v_0 only *)
regionFun = Function[{x, y}, Abs[vBipolar[x, y]] < v0];

(* Plot V/V_0 *)
Plot3D[
partialSum[5] @@ uvBipolar[x, y] // Re // Evaluate
, {x, -3, 3}, {y, -3, 3}
, Exclusions -> None
, RegionFunction -> regionFun
]


Next we take the gradient. Actually taking the derivative of a Fourier series is non-trivial. The $$u$$-derivative of $$\cos(nu)$$ introduces an extra factor of $$n$$, and if the coefficients do not go to zero fast enough, the term-by-term derivative will not converge. In this case though, the coefficients do go to zero fast enough for us to take a term-by-term derivative.

First define $$\partial V / {\partial u}$$ and $$\partial V / {\partial v}$$:

(* Derivatives of expansion terms *)
termUDerivative[n_][u_, v_] := D[term[n][u, v], u] // Evaluate;
termVDerivative[n_][u_, v_] := D[term[n][u, v], v] // Evaluate;

(* Derivatives partial sum *)
partialSumUDerivative[nMax_][u_, v_] := Sum[termUDerivative[n][u, v], {n, 0, nMax}];
partialSumVDerivative[nMax_][u_, v_] := Sum[termVDerivative[n][u, v], {n, 0, nMax}];


Since the scale factors for both coordinates are equal, the gradient is given by $$\nabla V = \frac{1}{h} \left( \frac{\partial V}{\partial u} \,\mathbf{a}_u + \frac{\partial V}{\partial v} \,\mathbf{a}_v \right):$$

(* Partial sum for gradient of V/V_0 *)
1 / hBipolar[u, v] * Plus[
partialSumUDerivative[nMax][u, v] uVectorBipolar[u, v],
partialSumVDerivative[nMax][u, v] vVectorBipolar[u, v]
];


Finally we can plot the gradient (I chose StreamPlot because VectorPlot arrows are too small):

Show[
ContourPlot[
partialSum[5] @@ uvBipolar[x, y] // Re // Evaluate
, {x, -3, 3}, {y, -3, 3}
, AspectRatio -> Automatic
, Exclusions -> None
, RegionFunction -> regionFun
],
StreamPlot[
gradientPartialSum[5] @@ uvBipolar[x, y] // Re // Evaluate
, {x, -3, 3}, {y, -3, 3}
, RegionFunction -> regionFun
]
]


• Thanks so much. If firstly we want to differentiate with respect to v from potential V and then plot potential and gradient function is possible to do it? If so the common differentiate should be applied?
– iman
Sep 12, 2020 at 12:59
• @iman I don't understand what you are asking. The StreamPlot is already a plot of the gradient. Sep 12, 2020 at 15:17
• Thanks. I have provided my question at the end of the body of the main primary question. Please see it and help again.
– iman
Sep 13, 2020 at 5:39
• @iman That is not a Mathematica problem, but a PDE problem (which is off-topic for this site). Firstly you might need different eigenfunctions and eigenvalues to satisfy $\partial V/ {\partial v} = 0$ at $v = v_0$. Secondly this is an homogeneous boundary condition. You will need a separate inhomogeneous boundary condition to determine the coefficients of the Fourier series. Sep 13, 2020 at 13:22
• Dear friend. I have been working on it, but I could not solve it. Could you please help me more in this regard? Thanks
– iman
Sep 17, 2020 at 21:15