4
$\begingroup$

I want to find Eigenvalues of following Sturm-Liouville problem

D[(2 - MA/κ[s]) q'[s], s] + W Ω^2 q[s] == 0,

Where Ω^2 is the eigenvalue I want to find.

To do so I use the following code

 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++,   
            Λ = Sqrt[(Cosh[L/(2 l)] - λ^2 + 
            ( λ^2 - 1) Cosh[s L/l])/(λ^2 ( Cosh[L/(2 l)] - 1))]; 
             X2[[j, 2]] = κ /. 
             FindRoot[χ (Λ^4) /κ^2 - χ + 
             2 Log[κ] + 2 L/(Pi 60000000) Cos[Pi s] == 0, {κ,
             1}]; X2[[j, 1]] = s; s = s + step;];
         Return[X2];]
 step1 = 1/10000;
 stepnumb = (1/2 + 1/2)/step1 + 1;


  SOL45 = Table[{0, 0}, {101}] ;
For [i = 1, i <= 101, i++,
 λ = 1.5;
 j = (i - 1)*0.524 + 1; X3 = Table[{0, 0}, {stepnumb}] ; 
  X3 = KS[1*60000000*Pi, λ, 1*60000000*Pi/6, j/1000, -1/2, 
    step1, stepnumb];
f = Interpolation[X3, Method -> "Spline", InterpolationOrder -> 20];
 ζ  = 3;
  Ct = 5/100; 
 χ = j/1000; 
 MA = N[(2 ζ  χ )  /(ζ + 1) * Ct ]; 
ρ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*60000000*Pi];
 l   = N[1*60000000*Pi/6]; 
 H = 60000000;   
Λ = 
Sqrt[(Cosh[
   L/(2 l)] - λ^2 + ( λ^2 - 1) Cosh[
     s L/l])/(λ^2 ( Cosh[L/(2 l)] - 1))]; κ[s] = 
f[s];
SOL45[[i, 1]] = j;
SOL45[[i, 2]] = 
NDEigenvalues[ {-D[(2 - MA/κ[s]) q'[s], s] /W, 
DirichletCondition[q[s] == 0, s == -1/2], 
DirichletCondition[ q[s] == 0, s == 1/2]}, q, {s, -1/2, 1/2}, 2];
]
      SOLPLOT45 = Table[{0, 0}, {100}];
      For[i = 1, i <= 100, i++,
       SOLPLOT45[[i, 1]] = SOL45[[i + 1, 1]]/1000;
 SOLPLOT45[[i, 2]] = Sqrt[SOL45[[i + 1, 2, 1]]];
  ]
 F45 = Interpolation[SOLPLOT45, InterpolationOrder -> 10]
P45 = Plot[F45[χ], {χ, 0, 0.0524}, 
 PlotStyle -> {Black, Dotted}, AspectRatio -> 2/3]

As a result I get the following set, whcih contains complex numbers and Plot wich doesn't make any sense, as I know that all Eigenvalues are real and positive in this case

enter image description here However the code works fine for other case, which just seems a bit different from expected Eigenvalue

 SOL1 = Table[{0, 0}, {101}] ;
For [i = 1, i <= 101, i++,
 λ = 1;
 j = (i - 1)*3 + 1; X3 = Table[{0, 0}, {stepnumb}] ; 
X3 = KS[1/4*60000000*Pi, λ, 1/4*60000000*Pi/6, j/1000, -1/2, 
 step1, stepnumb];
  f = Interpolation[X3, Method -> "Spline", InterpolationOrder -> 20];
 ζ  = 3; Ct = 2/100; χ = j/1000; 
 MA = N[(2 ζ  χ )  /(ζ + 1) * Ct ]; ρ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];
 l = N[1/4*60000000*Pi/6]; 
 H = 60000000;   Λ = 
 Sqrt[(Cosh[
     L/(2 l)] - λ^2 + ( λ^2 - 1) Cosh[
      s L/l])/(λ^2 ( Cosh[L/(2 l)] - 1))]; κ[s] = 
 f[s];
  SOL1[[i, 1]] = j;
 SOL1[[i, 2]] = 
 NDEigenvalues[ {-D[(2 - MA/κ[s]) q'[s], s] /W, 
  DirichletCondition[q[s] == 0, True]}, q, {s, -1/2, 1/2}, 2];
 ]
 SOLPLOT1 = Table[{0, 0}, {100}];
For[i = 1, i <= 100, i++,
SOLPLOT1[[i, 1]] = SOL1[[i + 1, 1]]/1000;
SOLPLOT1[[i, 2]] = Sqrt[SOL1[[i + 1, 2, 1]]];
]
 F1 = Interpolation[SOLPLOT1, InterpolationOrder -> 10]
 P1 = Plot[F1[χ], {χ, 0, 0.301}, PlotStyle -> {Black}, 
  AspectRatio -> 2/3]

enter image description here

Simplified code

For example if I take the following values

     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++,   
        Λ = 
       Sqrt[(Cosh[
         L/(2 l)] - λ^2 + ( λ^2 - 1) Cosh[
          s L/l])/(λ^2 ( Cosh[L/(2 l)] - 1))]; 
         X2[[j, 2]] = κ /. 
         FindRoot[χ (Λ^4) /κ^2 - χ + 
         2 Log[κ] + 2 L/(Pi 60000000) Cos[Pi s] == 0, {κ,
         1}]; X2[[j, 1]] = s; s = s + step;];
        Return[X2];]

   step1 = 1/10000;
  stepnumb = (1/2 + 1/2)/step1 + 1

 λ = 1;
 j = (79 - 1)*0.524 + 1; X3 = Table[{0, 0}, {stepnumb}] ; X3 = 
KS[1*60000000*Pi, λ, 1*60000000*Pi/6, j/1000, -1/2, step1, 
stepnumb];
f = Interpolation[X3, Method -> "Spline", InterpolationOrder -> 20];
ζ  = 3; Ct = 2/100; χ = j/1000; MA = 
N[(2 ζ  χ )  /(ζ + 1) * Ct ]; ρ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*60000000*Pi];
  l = N[1*60000000*Pi/6]; H = 60000000;   Λ = 
  Sqrt[(Cosh[
    L/(2 l)] - λ^2 + ( λ^2 - 1) Cosh[
     s L/l])/(λ^2 ( Cosh[L/(2 l)] - 1))]; κ[s] = f[s];

SOL100 = NDEigenvalues[ {-D[(2 - MA/κ[s]) q'[s], s] /W, 
 DirichletCondition[q[s] == 0, s == -1/2], 
  DirichletCondition[ q[s] == 0, s == 1/2]}, q, {s, -1/2, 1/2}, 2];
SOL100

Gives the answer {4.28055 - 15.1245 I, 4.28055 + 15.1245 I}

Where I know that the answer for the first Eigenvalue should be around 1.58

May I ask you for any ideas how to improve code to find Eigenvalues I do seek.

$\endgroup$
  • 1
    $\begingroup$ Readers typically avoid such long and slow-running blocks of code. Please consider simplifying your problem, probably by providing code to derive a single eigenvalue and comparing it with your expected value. $\endgroup$ – bbgodfrey Jul 23 '16 at 23:43
  • $\begingroup$ Followed your advice and added simplified case. $\endgroup$ – Alexander Jul 24 '16 at 9:39
  • $\begingroup$ I ran your simplified code using Version 10.4.1 for Microsoft Windows (64-bit) (April 11, 2016) and obtained {5.83485, 7.91617} as the answer. This is not the answer you predicted, but it is real. Please copy your simplified code into a new notebook, run it in a new session, and see what you get. $\endgroup$ – bbgodfrey Jul 24 '16 at 14:13
  • 1
    $\begingroup$ I recommend that you precede your code with Clear[Evaluate[Context[] <> "*"]] to prevent old values of constants from corrupting the output of KS. $\endgroup$ – bbgodfrey Jul 24 '16 at 16:00
9
$\begingroup$

For the "Simplified Code" provided in the OP's last edit, Λ in the definition of KS is identically 1, and the argument of FindRoot in the definition of KS is

frarg = χ (Λ^4) /κ^2 - χ + 2 Log[κ] + 2 L/(Pi 60000000) Cos[Pi s]
(* -0.041872 + 0.041872/κ^2 + 2. Cos[π s] + 2 Log[κ] *)

Zeros of frarg are shown in

ContourPlot[frarg == 0, {s, -1/2, 1/2}, {κ, 0, 1}]

enter image description here

and which solutionFindRoot chooses depends on the initial guess provided. For instance, an initial guess of 0.3 leads to the upper curve in the plot above, and corresponding eigenvalues {3.46382, 12.5554}, while an initial guess of 0.1 leads to the lower curve in the plot above, and corresponding eigenvalues {2.11177, 6.51226}. The latter is closer to the result predicted by the OP. On the other hand, the initial guess 1 used in the question leads to a mixture of solutions from the two branches. Consequently, the resulting eigenvalues are meaningless in that case.

For completeness, note that Solve also can be used to find κ. Since the two curves in the plot above do not cross, one need only choose the First or Last root from Solve to obtain the desired curve.

It seems likely that the erratic behavior in the first plot in the question occurs because of problems with κ.

Addendum: Replace KS, etc. by FunctionInterpolation

The "Simplified Code" section of the question can itself be simplified significantly by using FunctionInterpolation to replace KS, Interpolation and associated variables.

ClearAll[Evaluate[Context[] <> "*"]]
j = (79 - 1)*524/1000 + 1; χ = j/1000;
L = 60000000*Pi; l = 60000000*Pi/6; H = 60000000;  
λ = 1; ζ  = 3; Ct = 2/100; MA = (2 ζ χ)/(ζ + 1)*Ct; 
Λ = Sqrt[(Cosh[L/(2 l)] - λ^2 + (λ^2 - 1) Cosh[s L/l])/(λ^2 (Cosh[L/(2 l)] - 1))];
f = FunctionInterpolation[First[κ /. Quiet@NSolve[χ Λ^4/κ^2 - χ + 2 Log[κ] + 
    2 L/(Pi 60000000) Cos[Pi s] == 0, κ]], {s, -1/2, 1/2}, InterpolationPoints -> 101];
ρe = Exp[-L/(Pi H) Cos[Pi s]]; 
W = (Pi^2 Λ^4 )/(ζ + 1) (2 f[s]^2 ζ - (MA - 2) ρe )/(2 f[s] - MA); 
SOL100 = NDEigenvalues[ {-D[(2 - MA/f[s]) q'[s], s]/W , 
    DirichletCondition[q[s] == 0, s == -1/2], DirichletCondition[q[s] == 0, s == 1/2]}, 
    q, {s, -1/2, 1/2}, 2]
(* {2.11177, 6.51225} *)

(Use Last instead of First to choose the upper instead of the lower curve in the plot above.) Because f[s] is a slowly varying function, the number of evaluations of κ can be reduced from 10001 to 101, and InterpolationOrder reduced to its default value of 3.

Second Addendum: Curves for the original question

Looping over the code just above is straightforward and give the curves desired in the first part of the question. (Note that values of λ and Ct differ from those in the code just above so as to agree with the values in the first part of the question.

L = 60000000*Pi; l = 60000000*Pi/6; H = 60000000;  λ = 1.5; ζ  = 3; Ct = 5/100;
Λ = Sqrt[(Cosh[L/(2 l)] - λ^2 + (λ^2 - 1) Cosh[s L/l])/(λ^2 (Cosh[L/(2 l)] - 1))];
ρe = Exp[-L/(Pi H) Cos[Pi s]]; 
sol = Table[χ = j/1000; MA = (2 ζ χ) /(ζ + 1) * Ct; 
    f = FunctionInterpolation[First[κ /. 
    Quiet@NSolve[χ (Λ^4) /κ^2 - χ + 2 Log[κ] + 2 L/(Pi 60000000) Cos[Pi s] == 0, κ]], 
    {s, -1/2, 1/2}, InterpolationPoints -> 101];
    W = (Pi^2 Λ^4)/(ζ + 1) (2 f[s]^2 ζ - (MA - 2) ρe )/(2 f[s] - MA); 
    {χ, Sqrt @@ NDEigenvalues[{-D[(2 - MA/f[s]) q'[s], s]/W, 
    DirichletCondition[q[s] == 0, s == -1/2], DirichletCondition[q[s] == 0, s == 1/2]}, 
    q, {s, -1/2, 1/2}, 1]},
    {j, 1, 524/10, 524/1000}];
ListLinePlot[sol]

enter image description here

As noted in the first addendum, just replace First by Last to obtain results for the upper branch of κ in the first plot. Doing so yields

enter image description here

Finally, using L = 60000000*Pi/4; l = 60000000*Pi/24 along with λ = 1; Ct = 2/100 and a broader range of computation, {j, 1, 7492/25, 524/100}, reproduces the final plot in the question, validating the approach used in this answer. (Note that, for j much larger than 7492/25, no real values for κ exist near s == 0 for the parameters chosen, and FunctionInterpolation issues a plethora of error messages.)

$\endgroup$
  • $\begingroup$ Thank you for your comment I will try this method. $\endgroup$ – Alexander Jul 24 '16 at 17:54
6
$\begingroup$

I'll start with the equation the OP uses in the FindRoot[] call. It can be simplified and solved exactly symbolically. The following eq1, eq2, and eq3 are equivalent under the indicated transformations:

Clear[a, b, s, κ, L, χ, Λ];

cvt = {a -> 2 L/(Pi 60000000) Cos[Pi s] - χ,  (* convert to OP's "Simplified coode" *)
       b -> χ (Λ^4)};
eq1 = χ (Λ^4)/κ^2 - χ + 2 Log[κ] + 2 L/(Pi 60000000) Cos[Pi s] == 0;  (* original *)
eq2 = b/κ^2 + 2 Log[κ] + a == 0;                                      (* basic form *)
eq3 = b/κ2 + Log[κ2] + a == 0;                                        (* κ2 == κ^2 *)

Simplify[eq1] ===
 Simplify[eq2 /. cvt] ===
 Simplify[eq3 /. κ2 -> κ^2 /. cvt, κ > 0]
(*  True  *)

All of the Solve[] commands below give the warning:

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>

The equation eq2, like the OP's eq1, has two solutions, but clearly both cannot be correct, since they have opposite signs and the equation is not symmetric with respect to κ:

Solve[b/κ^2 + 2 Log[κ] + a == 0, κ]
(*
  {{κ -> -((I Sqrt[b])/Sqrt[ProductLog[-b E^a]])},
   {κ -> (I Sqrt[b])/Sqrt[ProductLog[-b E^a]]}}
*)

If we replace κ^2 by κ2 (eq3), we get just one solution.

Solve[b/κ2 + Log[κ2] + a == 0, κ2]    
(*  {{κ2 -> -(b/ProductLog[-b E^a])}}  *)

This is the solution we'll start with. The warning by Solve[] suggests trying other branches of ProductLog[], that is, ProductLog[r, -b E^a] for various integers r. To verify the solutions for the OP's equation, we will need to plug in the actual numerical parameters for the system. This will be done below. It turns out there are two branches that solve the equation, r == 0 and r == -1. We use this code to construct the function κ[s, r] in the solution further down.

ksol = κ -> Sqrt[κ2] /. Solve[b/κ2 + Log[κ2] + a == 0, κ2] /. cvt
(*
  {κ -> Sqrt[-((Λ^4 χ)/ProductLog[-E^(-χ + (L Cos[π s])/(30000000 π)) Λ^4 χ])]}
*)

Here is the adapted code. I've made some changes to make the interdependencies a bit easier to manage.

(* initialize numeric parameters *)
step1 = 1/10000;
stepnumb = (1/2 + 1/2)/step1 + 1;

λ = 1;
j = (79 - 1)*0.524 + 1;
ζ = 3;
Ct = 2/100;
χ = j/1000;
MA = N[(2 ζ χ)/(ζ + 1)*Ct];
L = 1*60000000*Pi;
l = 1*60000000*Pi/6;
H = 60000000;

(* initialize functions *)
ClearAll[X2, X3, SOL100, κ, W, Λ, ρe];
KS[L_, λ_, l_, χ_, s0_, step0_, stepnumb_, root_] := 
  Module[{κ, Λ},
   Table[
    With[{s = s0 + (j - 1) step0},
     Λ = Sqrt[(Cosh[L/(2 l)] - λ^2 + (λ^2 - 1) Cosh[s L/l])/(λ^2 (Cosh[L/(2 l)] - 1))];
     {s,  
      κ /. FindRoot[χ (Λ^4)/κ^2 - χ + 2 Log[κ] + 2 L/(Pi 60000000) Cos[Pi s] == 0, 
        Evaluate@If[root == 0, {κ, 1, 0.35}, {κ, 0.1, 0.2}]]}
     ],
    {j, stepnumb}]
   ];

mem : X2[root_] := 
  mem = KS[1*60000000*Pi, λ, 1*60000000*Pi/6, j/1000, -1/2, step1, stepnumb, root];

ρe[s_] := Exp[-L/(Pi H) Cos[Pi s]];
W[s_, r_] := (Pi^2 Λ[s]^4)/(ζ + 1) (2 κ[s, r]^2 ζ - (MA - 2) ρe[s])/(2 κ[s, r] - MA);
Λ[s_] := Sqrt[(Cosh[L/(2 l)] - λ^2 + (λ^2 - 1) Cosh[s L/l])/(λ^2 (Cosh[L/(2 l)] - 1))];
κ[s_, r_] := 
  Sqrt[-((Λ[s]^4 χ)/ProductLog[r, -E^(-χ + (L Cos[π s])/(30000000 π)) Λ[s]^4 χ])];

mem : SOL100[root_] := 
  mem = NDEigenvalues[{-D[(2 - MA/κ[s, root]) q'[s], s]/W[s, root],
     DirichletCondition[q[s] == 0, s == -1/2],
     DirichletCondition[q[s] == 0, s == 1/2]},
    q, {s, -1/2, 1/2}, 2];

The two solutions, root = 0 and root = -1:

Block[{root = 0},
 Print["eigenvalues"[root] -> SOL100[root]];
 Show[
  Plot[κ[s, root], {s, -1/2, 1/2}, 
   PlotStyle -> {Black, Thickness[0.015]},
   AxesOrigin -> {0, 0}],
  ListPlot[X2[root], PlotStyle -> Red]
  ]
 ]

(*  "eigenvalues"[0] -> {3.46382,12.5554}  *)

Mathematica graphics

Block[{root = -1},
 Print["eigenvalues"[root] -> SOL100[root]];
 Show[
  Plot[κ[s, root], {s, -1/2, 1/2}, 
   PlotStyle -> {Black, Thickness[0.015]},
   AxesOrigin -> {0, 0}],
  ListPlot[X2[root], PlotStyle -> Red]
  ]
 ]

(*  "eigenvalues"[-1] -> {2.11177,6.51226}  *)

Mathematica graphics

It turns out to be difficult to test the correctness of the solutions symbolically, so I settled for checking numerically. For values of r in ProductLog[r, z] other than 0 or -1, we get nonzero imaginary residuals that grow with Abs[r]. That indicated to me that there were no other solutions. Below we give code for checking the two valid cases. (Note ProductLog[0, z] === ProductLog[z].)

SetPrecision[
 eq1 /. Equal -> Subtract /. ksol /. Λ -> Λ[s],
 $MachinePrecision] /. 
   s -> RandomReal[{-1/2, 1/2}, 1000, WorkingPrecision -> $MachinePrecision] // Abs // Max

(*  0.*10^-15  *)

SetPrecision[
 eq1 /. Equal -> Subtract /. ksol /. {ProductLog[z_] :> ProductLog[-1, z]} /. Λ -> Λ[s],
 $MachinePrecision] /. 
   s -> RandomReal[{-1/2, 1/2}, 1000, WorkingPrecision -> $MachinePrecision] // Abs // Max

(*  0.*10^-14  *)
$\endgroup$
  • 1
    $\begingroup$ It should be noted that only two branches of the Lambert function have real values; $W_0(x)=W(x)$ is real for $x\ge-\frac1e$, and $W_{-1}(x)$ is real for is real for $-\frac1e\le x<0$. $\endgroup$ – J. M. will be back soon Jul 27 '16 at 14:41

Not the answer you're looking for? Browse other questions tagged or ask your own question.