3
$\begingroup$

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!

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Bongy May 11 '18 at 16:10
  • $\begingroup$ 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. $\endgroup$ – Bongy May 11 '18 at 16:19
  • $\begingroup$ 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}}] $\endgroup$ – Lotus May 14 '18 at 8:41
  • $\begingroup$ By the way if your main question is answered could you please accept my answer ? Thanks! $\endgroup$ – Lotus May 14 '18 at 8:42
2
$\begingroup$

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

enter image description here

Previously found minimum

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

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

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

enter image description here

$\endgroup$
1
$\begingroup$

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] == 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}]

enter image description here

$\endgroup$

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.