2

I am trying to find eigenvalues of a Sturm-Liouville equation using Mathematica10.2, but setting different precisions in ParametricNDSolve gives me different results for the eigenvalues. I start by solving a transcendental equation to get some parameters for my function

KS[L0_, λ0_, l0_, χ0_, s0_, step0_, stepnumb0_] := 
  Module[
      {L = L0, l = l0, λ = λ0, χ = χ0, j, s = s0, 
       step = step0, stepnumb = stepnumb0}, 
    X2 = Table[{0, 0}, {stepnumb}];
    For[j = 1, j <= stepnumb, j++,   
      Λ = 
        N[Sqrt[(Cosh[L/(2 l)] - λ^2 + ( λ^2 - 1) Cosh[s L/l])/
          (λ^2 ( Cosh[L/(2 l)] - 1))], 45]; 
      X2[[j, 2]] = 
        κ /. FindRoot[
               χ (Λ^4) /κ^2 - χ + 2 Log[κ] + 2 L/(Pi 60000000) Cos[Pi s] == 0, 
               {κ, 1}, 
               WorkingPrecision -> 40]; 
      X2[[j, 1]] = s; 
      s = s + step;];
    Return[X2];]
step1=1/10000
stepnumb = (1/2 + 1/2)/step1 + 1
X3 = Table[{0, 0}, {stepnumb}]; 
X3 = 
  KS[1/4*60000000*Pi, 1, 1/4*60000000*Pi/6, 200/1000, -1/2, step1, stepnumb];

ListPlot[X3]

Then I interpolate it.

f = Interpolation[X3, Method -> "Spline", InterpolationOrder -> 20];

I define different values

 ζ = 3; 
 Ct = 2/100; 
 χ = 200/1000; 
 MA = N[(2 ζ  χ )  /(ζ + 1) * Ct , 45]; 
 ρe = 
   Exp[-L/(Pi H) Cos[Pi s]]; W = (Pi^2 Λ^4 )/(ζ + 1) (2 κ[s]^2 ζ - (MA - 2) ρe )/(2 κ[s] - MA); 
 L = N[1/4*60000000*Pi, 45];  
 λ = 1; 
 l = N[1/4*60000000*Pi/6, 45]; 
 H = 60000000;
 Λ = 
   Sqrt[(Cosh[L/(2 l)] - λ^2 + ( λ^2 - 1) Cosh[s L/l])/(λ^2 (Cosh[L/(2 l)] - 1))]; 
 κ[s] = f[s];

Then I try to find the eigenvalues.

sol3 = 
  ParametricNDSolve[
    {-D[(2 - MA/κ[s]) q'[s], s] + W Ω ^2 q[s] == 0, q[-1/2] == 0, q[1/2] == 0}, 
    q, {s, -1/2, 1/2}, {Ω}, 
    Method -> {"Shooting"},  
    WorkingPrecision -> 35]

FindRoot[q[Ω][0.0] /. sol3, {Ω, 1 + 5/10, 1 + 7/10}, Method -> {"Secant"}]

But in return it gives completly random answers.

May I ask you for advice on finding real eigenvalues? For example how to find first 2-3 Eigenvalues.

Alexander
  • 77
  • 5

1 Answers1

1

The problem here is the method you are using for Findroot[]

Using Newton's method reliably gives the answer:

FindRoot[q[Ω][0.0] /. sol3, {Ω, 1 + 5/10}]

`{Ω -> 1.57047}

Note that with two starting values, you will always use the Secant method, and specifying so as you did was redundant.

Running the Findroot[] multiple times, I've noticed that sometimes, it yields an answer of ~1.39, though as you know the answer is $1.5\le \Omega \le 1.7$, you can ignore these.

Feyre
  • 8,597
  • 2
  • 27
  • 46
  • May I ask you if there is a way of finding for example first 3 eigenvalues of the problem? – Alexander Jul 21 '16 at 13:23
  • @Alexander Normally I'd use NSolve[], but I don't think that's applicable here – Feyre Jul 21 '16 at 13:28
  • Thank you for your help but unfortunatelly it still gives me some random answers if I do run it multiple times – Alexander Jul 21 '16 at 13:32
  • @Alexander It constantly gives me 1.57047, and sometimes 1.398 which I suspect is another Eigen value, regardless of the domain you've chosen. I've run the code around 20 times and those are the only numbers I get. Have you tried quitting the kernel and rerunning everything? – Feyre Jul 21 '16 at 13:38
  • Yes I did, and answer 1.57047, does not apper for me at all, but in fact I know that this answer should be correct one. – Alexander Jul 21 '16 at 13:40
  • 1
    maybe I should use commands like NDEigenvalues, but I am not sure if they are appropriate ones – Alexander Jul 21 '16 at 13:42
  • @Alexander Curious and curiouser. Then this is likely to be a version issue, try updating to 10.4 if you haven't already. – Feyre Jul 21 '16 at 13:44
  • Tried to upgrade, it gives mostly correct answers, however some of them are completely of the grid, still it is much better then it was before. – Alexander Jul 21 '16 at 16:47
  • 1
    You might try FindAllCrossings, if the roots sought are real. For complex roots, divide them out of q[Ω][0.0] /. sol3 as you find them, and they will not be found again. – bbgodfrey Jul 22 '16 at 14:38