0
$\begingroup$

I am trying to solve the following problem numerically and I faced two things: 1- It takes too long time to give solutions. 2- there are many roots for the equation but NSolve doesn't give all thee possible roots, it seems that it choose the first value that satisfy the equation! The solutions should be as in the fig below

λ11 = (-k^2 + Sqrt[x^2 - ub1^2])^(1/2);
λ12 = (-k^2 - Sqrt[x^2 - ub1^2])^(1/2);
f11 = (I λ11 - k)^2/(x - ub1);
f12 = (I λ12 - k)^2/(x - ub1);
λ21 = (-k^2 + Sqrt[x^2 - ub2^2])^(1/2);
λ22 = (-k^2 - Sqrt[x^2 - ub2^2])^(1/2);
g21 = (I λ21 + k)^2/(x + ub2);
g22 = (I λ22 + k)^2/(x + ub2);

M = {{1, 1, -1, -1}, {f11, f12, -g21, -g22}, {λ11, λ12, λ21, λ22}, 
    {λ11 f11, λ12 f12, λ21 g21, λ22 g22}};    
ub1:=0.25;
ub2:=0.25;
Table[NSolve[{Det[M]==0},x,WorkingPrecision->4],{k,-1,1,0.1}]

enter image description here

$\endgroup$
  • 4
    $\begingroup$ The speed issue is resolved by computing Det[M]outside Table: For instance, s = Numerator[Det[M] // FullSimplify] and x /. Table[NSolve[s == 0, x], {k, -1, 1, 0.1}]. However, it may be necessary to use FindRoot to find all the roots. $\endgroup$ – bbgodfrey Apr 2 '16 at 1:11
  • $\begingroup$ @bb, the resulting equations will be algebraic, so I think NSolve[] can handle it (though preprocessing with GroebnerBasis[] might be needed). $\endgroup$ – J. M. will be back soon Apr 2 '16 at 1:55
  • 2
    $\begingroup$ @J.M. Following you suggestion, as I understand it, I tried First@GroebnerBasis[Det[M], x] and was able to reproduce the curves in the question. Unfortunately, the solutions do not for the most part satisfy Det[M]==0. This and other attempts to obtain the curves lead me to believe that the equation in the question do not, in fact, yield these curves. $\endgroup$ – bbgodfrey Apr 2 '16 at 5:15
  • 2
    $\begingroup$ Faster if given exact input, that is, use ub1 = 1/4; ub2 = 1/4; and Table[NSolve[{Det[M] == 0}, x], {k, -1, 1, 1/10}]. Not sure what is the bottleneck offhand, might have to do with removing the parasite solutions that arise from having the variables in radicals. $\endgroup$ – Daniel Lichtblau Apr 2 '16 at 16:25
  • 1
    $\begingroup$ As requested, I added the GroebnerBasis plot. Please add to your question the code that you used to create your plot. Thanks. $\endgroup$ – bbgodfrey Apr 2 '16 at 20:46
2
$\begingroup$

As requested, here is plot obtained using GroebnerBasis (suggested by J.M.).

ub1 = 1/4; ub2 = 1/4;
λ11 = (-k^2 + Sqrt[x^2 - ub1^2])^(1/2);
λ12 = (-k^2 - Sqrt[x^2 - ub1^2])^(1/2);
f11 = (I λ11 - k)^2/(x - ub1);
f12 = (I λ12 - k)^2/(x - ub1);
λ21 = (-k^2 + Sqrt[x^2 - ub2^2])^(1/2);
λ22 = (-k^2 - Sqrt[x^2 - ub2^2])^(1/2);
g21 = (I λ21 + k)^2/(x + ub2);
g22 = (I λ22 + k)^2/(x + ub2);
M = {{1, 1, -1, -1}, {f11, f12, -g21, -g22}, {λ11, λ12, λ21, λ22}, 
     {λ11 f11, λ12 f12, λ21 g21, λ22 g22}};

s = FullSimplify[Det[M]];
GroebnerBasis[s, x][[1]]
(* 1 - 256 k^4 + 16384 k^8 - 128 x^2 - 18432 k^4 x^2 - 262144 k^8 x^2 + 6144 x^4 + 
   393216 k^4 x^4 + 1048576 k^8 x^4 - 131072 x^6 - 2097152 k^4 x^6 + 1048576 x^8 *)
x /. Solve[% == 0, x] // Flatten;
Plot[Evaluate[%], {k, -1, 1}, AspectRatio -> 1, Frame -> True, FrameLabel -> {k, x}]

enter image description here

It differs from the figure in the question (apart from colors) by the presence short curve segments in the corners. They can be eliminated, if desired, by using Evaluate[%[[5;;8]]] instead of Evaluate[%] in the final line of code.

$\endgroup$
  • $\begingroup$ These would be most likely from "parasite" solutions due to presence of radicals (probably you knew that but maybe it's not obvious to everyone). $\endgroup$ – Daniel Lichtblau Apr 2 '16 at 20:53
  • $\begingroup$ @bb, In fact, I eliminated those curves at the corners since I was interested to find solution for Abs[x]<0.4 and I forgot to add them to the plot in the question. I think we did the same thing, but I am interested to find these solutions using NSolve not Solve. As I mentioned earlier I would like to compare the results for Solve and NSolve $\endgroup$ – HD2006 Apr 2 '16 at 22:24

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.