# Numerically evaluating and plotting functions involving complete Elliptic integrals

I have an equation like this

$$-K(y)-\frac{\pi J}{2 \sqrt{\lambda }}+\frac{E(y)}{1-y}=0$$

1/(1 - y) EllipticE[y] - EllipticK[y] - (J π /(2 Sqrt[λ])) == 0


Where $E(y)$ and $K(y)$ are the Elliptic integrals in the Mathematica notation and $y$ is Real. I want to solve the equation numerically for $y$ for different values of the variables $(\lambda, J)$. That is $(\lambda, J)$ acts as a $(x,y)$ grid where $y$ gives the values along $z$ axis. Now I want to take these values of $y$ for different values of $(\lambda, J)$ and evaluate the function

$$\Delta= \frac{2 \sqrt{\lambda} }{\pi} \left(K(y)+\frac{\left(\sqrt{y}-1\right) E(y)}{1-y}\right)$$

(2 Sqrt[λ])/π ((Sqrt[y] - 1)/(1 - y) EllipticE[y] +
EllipticK[y])


and plot it in 3d along the $\Delta, \lambda, J$ axes. So the values on the $(\lambda, J)$ plane remains the same, but instead of $y$, we want to plot $\Delta(y)$ at each $(\lambda, J)$ value.

Can someone help me write the code? Thanks!

There were some syntax errors in your code, but here is a cleaned-up version:

Define the function (I used lam instead of Lambda):

fun[lam_, j_, y_] :=
1/(1 - y) EllipticE[y] - EllipticK[y] - (j Pi/(2 Sqrt[lam]))


Note that FindRoot gives the root - for example:

FindRoot[fun[1, 1, y] == 0, {y, 0.1}]

(* y -> 0.642898 *)


Use Table to generate results for different values of lam, j and y:

data = Flatten[
Table[{lam, j, y /. FindRoot[fun[lam, j, y] == 0, {y, 0.1}]}, {lam,
0.1, 0.5, 0.1}, {j, 0.1, 0.5, 0.1}], 1]


Define the Del function:

del[lam_, j_,
y_] := (2 Sqrt[lam])/
Pi ((Sqrt[y] - 1)/(1 - y) EllipticE[y] + EllipticK[y])


Take y from the previously generated results and generate a new data set:

deldata =
Flatten[Table[{lam, j, del[lam, j, y]}, {lam, 0.1, 0.5, 0.1}, {j,
0.1, 0.5, 0.1}, {y, data[[All, 3]]}], 2]


Use ListPointPlot3D to plot

ListPointPlot3D[deldata]


Please note that the values of lam and j I have used are completely fictitious and may not be what you would want to use at all. But I do hope that the procedure is clear.

• This works fine, thanks. One small problem, in the ListPlot the x,y axes are okay, but along the z axis the plot keeps showing values of y instead of values of \Del(y). Any idea why? – Bongy May 11 '18 at 16:10
• What I mean to say that ListPointPlot3D[deldata] and ListPlot3D[deldata] have different values on the z axis. One shows values of \Del(y) and other shows values of y. I'm bit confused there. – Bongy May 11 '18 at 16:19
• I think because of the values that I chose the roots are quite close to each other and therefore ListPlot3D doesn't really show this difference. Compare ListPointPlot3D[{{0, 0, 1}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1.01}, {1.01, 0, 0}, {0, 1.01, 0}}] and ListPlot3D[{{0, 0, 1}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1.01}, {1.01, 0, 0}, {0, 1.01, 0}}] – Lotus May 14 '18 at 8:41
• By the way if your main question is answered could you please accept my answer ? Thanks! – Lotus May 14 '18 at 8:42

Here another approach.

Equation 1 depends only on J/Sqrt[lam] , I call that Jlam.

eq1 = -((J \[Pi])/(2 Sqrt[lam])) + EllipticE[y]/(1 - y) -
EllipticK[y] == 0 /. J/Sqrt[lam] -> Jlam

eq2 = delta == (2 Sqrt[lam])/\[Pi] ((Sqrt[y] - 1)/(1 - y) EllipticE[
y] + EllipticK[y])


Get an impression of solution with

ContourPlot[Evaluate[eq1], {Jlam, -3, 3}, {y, -4, 4},
WorkingPrecision -> 20]


Show, how Jlam depends on y

sol1 = First@Solve[eq1, Jlam]

Plot[Jlam /. sol1, {y, -30, 1}, PlotRange -> {-1, 1},
Epilog -> {Red, Point[{yMin, JlamMin}]}] Previously found minimum

mi = Minimize[{Jlam /. sol1, -10 < y < 0}, y] // Quiet // FullSimplify

{JlamMin, yMin} = {mi[], y /. mi[]}


Get y as function of Jlam with FindRoot. Here I regard only the right wing of Jlam near zero.

 yfr[Jlam_?NumericQ] :=
y /. First@
FindRoot[-((Jlam \[Pi])/2) + EllipticE[y]/(1 - y) - EllipticK[y] ==
0, {y, -.1, yMin, .99}]


Now calculate delta with y from yfr

delta[J_?NumericQ, lam_?NumericQ] := {y = yfr[J/Sqrt[lam]], (
2 Sqrt[lam] (((-1 + Sqrt[y]) EllipticE[y])/(1 - y) +
EllipticK[y]))/\[Pi]}

Plot3D[Last@delta[J, lam], {J, 0, 3}, {lam, 0, 4},
ImageSize -> 300] // Timing Here's an alternate method to produce @Akku's answer. Again, the dependence on $J$ and $\sqrt{\lambda}$ can be replaced by a dependence on a single parameter $x$:

eqn = 1/(1-y) EllipticE[y] - EllipticK[y] - (J \[Pi]/(2 Sqrt[λ])) == 0;
eqn = eqn /. {y->y[x], J -> x Sqrt[λ]};

eqn //TeXForm


$-K(y(x))+\frac{E(y(x))}{1-y(x)}-\frac{\pi x}{2}=0$

Next, we can turn this equation into an ODE:

ode = D[eqn, x];

ode //TeXForm


$\frac{y'(x) (E(y(x))-K(y(x)))}{2 (1-y(x)) y(x)}-\frac{y'(x) (E(y(x))-(1-y(x)) K(y(x)))}{2 (1-y(x)) y(x)}+\frac{y'(x) E(y(x))}{(1-y(x))^2}-\frac{\pi }{2}=0$

We can solve this ODE using NDSolve, but we will need an initial condition. Note that the equation evaluated at $y(x)=0$ gives:

eqn /. y[x]->0


-(π x)/2 == 0

So, the initial condition is $y(0)=0$. Now, we are ready to use NDSolveValue:

sol = NDSolveValue[{ode, y == 0}, y, {x, -1, 1000}];


NDSolveValue::ndsz: At x == -0.308082, step size is effectively zero; singularity or stiff system suspected.

We can use the above solution to create a function $yy(J, \lambda)$:

yy[J_, λ_] := sol[J/Sqrt[λ]]


Finally, we can insert this function into your $\Delta$ equation:

Δ = (2 Sqrt[λ])/π ((Sqrt[yy[J,λ]] - 1)/(1 - yy[J,λ]) EllipticE[yy[J,λ]] + EllipticK[yy[J,λ]]);


Visualization:

Plot3D[Δ, {J,0,3}, {λ,0,4}] 