2
$\begingroup$

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]

$\endgroup$
7
  • $\begingroup$ Could you check your second definition (V)? It is not written properly. $\endgroup$
    – Sumit
    Sep 10 '20 at 11:41
  • $\begingroup$ Thanks. I think it is true. Anyway, I provide again it in the original form as a picture. $\endgroup$
    – iman
    Sep 10 '20 at 12:04
  • $\begingroup$ 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. $\endgroup$
    – Sumit
    Sep 10 '20 at 13:31
  • $\begingroup$ 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 $\endgroup$
    – iman
    Sep 10 '20 at 13:48
  • $\begingroup$ 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. $\endgroup$
    – iman
    Sep 10 '20 at 18:38
2
$\begingroup$

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

Fourier series for boundary condition

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
]

Plot of function

Gradient

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 *)
gradientPartialSum[nMax_][u_, v_] :=
  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
    , ContourShading -> None
    , Exclusions -> None
    , RegionFunction -> regionFun
  ],
  StreamPlot[
    gradientPartialSum[5] @@ uvBipolar[x, y] // Re // Evaluate
    , {x, -3, 3}, {y, -3, 3}
    , RegionFunction -> regionFun
  ]
]

Plot of gradient

$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – iman
    Sep 12 '20 at 12:59
  • $\begingroup$ @iman I don't understand what you are asking. The StreamPlot is already a plot of the gradient. $\endgroup$
    – yawnoc
    Sep 12 '20 at 15:17
  • $\begingroup$ Thanks. I have provided my question at the end of the body of the main primary question. Please see it and help again. $\endgroup$
    – iman
    Sep 13 '20 at 5:39
  • $\begingroup$ @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. $\endgroup$
    – yawnoc
    Sep 13 '20 at 13:22
  • $\begingroup$ Dear friend. I have been working on it, but I could not solve it. Could you please help me more in this regard? Thanks $\endgroup$
    – iman
    Sep 17 '20 at 21:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.