I wish to solve an equation for the field lines of an electric field. This is the textbook question 1.3.8(a) from D. Dubin’s Numerical and Analytical Methods for Scientists and Engineers, Using Mathematica (ISBN-13: 978-0471266105), however it is not an assigned homework problem (so feel free to give direct answers!)

Theory

Given the potential ϕ(ρ,z)=z/(ρ2+z2)3/2ϕ(ρ,z)=z/(ρ2+z2)3/2\phi(\rho,z)=z/(\rho^2+z^2)^{3/2}, which has a corresponding electric field

E(ρ,z)=3ρz(ρ2+z2)5/2ˆρ+(3z2(ρ2+z2)5/2−1(ρ2+z2)3/2)ˆz\mathbf{E}(\rho,z) = {3\rho z \over(\rho^2 + z^2)^{5/2}}\mathbf{\hat{\rho}} + \left({3z^2 \over(\rho^2 + z^2)^{5/2}}-{1\over(\rho^2 + z^2)^{3/2}}\right) \mathbf{\hat{z}}

as given by E=−∇ϕ\mathbf{E}=-\mathbf{\nabla}\phi, the equation for a field line at position ss is given by

dρds=Eρ(ρ,z)|Eρ(ρ,z)|.{\mathrm{d}\rho\over\mathrm{d}s}={E_\rho(\rho,z)\over|E_\rho(\rho,z)|}.

This is defined in my code as ODEÏ[s], with the magnitude |Eρ(ρ,z)||E_\rho(\rho,z)| defined as Emag[Ï[s], z[s]].

Mathematica code

Here I am tryin to solve the equation for the ρ\rho-component only.

ClearAll[“Global`*”];

ODEÏ[s_] =

Ï'[s] == (3 Ï[s] z[s])/(Ï[s]^2 + z[s]^2)^(5/2) 1/Emag[Ï[s], z[s]];

Emag[Ï[s_], z[s_]]

= Sqrt[((3 Ï[s] z[s])/(Ï[s]^2 + z[s]^2)^(5/2))^2

+ ((3 z[s]^2)/(Ï[s]^2 + z[s]^2)^(5/2) – 1/(Ï[s]^2 + z[s]^2)^(3/2))^2];

DSolve[ODEÏ[s], Ï[s], s]

which produces the following error message:

Solve::ifun: Inverse functions are being used by Solve, so some

solutions may not be found; use Reduce for complete solution

information. >>

then replicates the DSolve statement without solving anything.

Could you please include in your response what DSolve is trying to do to solve the equation and why the error has been produced, perhaps with a quick description of why your solution works.

=================

=================

1 Answer

1

=================

The simple answer (Dec. 12, 15)

After a rather long journey with Mathematica I found out that the result can be simplified drastically to fairly simple textbook wisdom.

Introducing polar coordinates {r,a} in the {ρ\rho, z} plane, which means in fact working in spherical coordinates, we have the following results

Lines of equal potential

ϕ=z(ρ2+z2)3/2=cos(a)r2=const\phi =\frac{z}{\left(\rho ^2+z^2\right)^{3/2}}=\frac{\cos (a)}{r^2}=\text{const}

r(a)=c√cos(a)r(a)=c \sqrt{\cos (a)}

Field lines

ψ=z2(ρ2+z2)3/2=sin2(a)r=const\psi =\frac{z^2}{\left(\rho ^2+z^2\right)^{3/2}}=\frac{\sin ^2(a)}{r}=\text{const}

r(a)=csin2(a)r(a)=c \sin ^2(a)

Overview with simple streamplot

As a first step, disregarding the cylindrical geometry, you could use StreamPlot:

Ï•[Ï_, z_] = z/(Ï^2 + z^2)^(3/2);

fE = {-D[Ï•[Ï, z], Ï], -D[Ï•[Ï, z], z]}

(*

{(3 z Ï)/(z^2 + Ï^2)^(5/2), (3 z^2)/(z^2 + Ï^2)^(5/2) – 1/(z^2 + Ï^2)^(3/2)}

*)

StreamPlot[fE, {Ï, 0.1, 4}, {z, -4, 4}]

Calculation of the field lines

In cylindrical coordinates with azimuthal symmtery the field lines ρ(z)\rho (z) are determined by the ODE

dρ/Eρ=dz/Ezd\rho/ E _\rho = dz / E _z

Potantial and field strength are, respectively,

Ï•[Ï_, z_] = z/(Ï^2 + z^2)^(3/2);

fE = -{D[Ï•[Ï, z], Ï], D[Ï•[Ï, z], z]} //

Simplify

(*

Out[3]= {(3 z Ï)/(z^2 + Ï^2)^(5/2), (

2 z^2 – Ï^2)/(z^2 + Ï^2)^(5/2)}

*)

With the field vector fE we have the ODE

eq = Ï'[z] == (fE[[1]]/fE[[2]] /. Ï -> Ï[z])

(*

Out[5]= Derivative[1][Ï][z] == (3 z Ï[z])/(2 z^2 – Ï[z]^2)

*)

The solution is

Ïs0 = Ï[z] /. DSolve[eq, Ï[z], z];

Length[Ïs0]

(* Out[148]= 6 *)

It consists of 6 componentes and containes one Integration constant C1.

At z = 0 it becomes

Simplify[Ïs0 /. z -> 0, C[1] âˆˆ Reals]

(*

Out[149]= {-E^C[1], E^C[1], 0, 0, 0, 0}

*)

Hence we have identified C1 -> Log[Ï0] and the solution becomes

Ïs[z_, Ï0_] = Ïs0 /. C[1] -> Log[Ï0];

An appropriate subset of these 6 solutions is the explicit equation for the field lines.

As an example the first component is

Ïs[z, Ï0][[1]];

−√ρ023−z2+3√2ρ06+27ρ02z4−18ρ04z2+3√3√27ρ04z8−4ρ06z633√2−23√2ρ02z23√2ρ06+27ρ02z4−18ρ04z2+3√3√27ρ04z8−4ρ06z6+3√2ρ0433√2ρ06+27ρ02z4−18ρ04z2+3√3√27ρ04z8−4ρ06z6-\sqrt{\frac{\text{$\rho $0}^2}{3}-z^2+\frac{\sqrt[3]{2 \text{$\rho $0}^6+27 \text{$\rho $0}^2 z^4-18 \text{$\rho $0}^4 z^2+3 \sqrt{3} \sqrt{27 \text{$\rho $0}^4 z^8-4 \text{$\rho $0}^6 z^6}}}{3 \sqrt[3]{2}}-\frac{2 \sqrt[3]{2} \text{$\rho $0}^2 z^2}{\sqrt[3]{2 \text{$\rho $0}^6+27 \text{$\rho $0}^2 z^4-18 \text{$\rho $0}^4 z^2+3 \sqrt{3} \sqrt{27 \text{$\rho $0}^4 z^8-4 \text{$\rho $0}^6 z^6}}}+\frac{\sqrt[3]{2} \text{$\rho $0}^4}{3 \sqrt[3]{2 \text{$\rho $0}^6+27 \text{$\rho $0}^2 z^4-18 \text{$\rho $0}^4 z^2+3 \sqrt{3} \sqrt{27 \text{$\rho $0}^4 z^8-4 \text{$\rho $0}^6 z^6}}}}

In order for the square root appearing in the solution to stay real we must have

cond = – ((2 Ï0)/(3 Sqrt[3])) < z < (2 Ï0)/(3 Sqrt[3]);
Now we are ready to visualize the result:
A plot of the field lines is then accomplished by this function
pE[Ï0_] :=
Plot[Re[Take[Ïs[z, Ï0], {1, 4}]], {z, -((2 Ï0)/(3 Sqrt[3])), (
2 Ï0)/(3 Sqrt[3])}, PlotStyle -> {Thin, Red}, PlotRange -> All,

PlotLabel -> “Field lines (red) and equiotentials (colored)”,

AxesLabel -> {“z”, “Ï”}]

pE[2]

A plot of the equipotetial lines is

Ï•[Ï_, z_] = z/(Ï^2 + z^2)^(3/2);

pÏ• = ContourPlot[Ï•[Ï, z], {z, -5, 5}, {Ï, -4, 4},

Contours -> 15, Frame -> True , Axes -> True, AxesLabel -> {“z”, “Ï”}]

A joint plot shows both line families

Show[{pÏ•, pE[#] & /@ {1, 2, 3, 4, 8, 16}},

PlotLabel -> “Field lines (red) and equipotentials (colored)”]

Although this is certainly a useful technique for producing/viewing field lines, the problem I am hoping to solve has more to with finding a solution to the ODE rather than producing valid output through other techniques. It’s more of a math excercise than a computing one, I guess.

– Harry Smith

Dec 9 ’15 at 22:49

@Harry Smith: in the meantime I have also solved the ODE and found the explicit formula for the field lines. I shall post it soon. I consider this a MMA problem, and do not agree with closing the question.

– Dr. Wolfgang Hintze

Dec 11 ’15 at 12:44