Stationary solution of PDE and eigenvalues

Question

I need to find the eigenvalues of the stationary solution ($\dot{u}=0$) of the following PDE

$ \dot{u}=-d u + (d-2)x u'-v(d)L_{d}^{1}(u'+2x u'') $ with $ L_{d}^{n}(f)=\int_0^\infty dy~y^{-1+d/2}\left(\frac{-2 a(1+y)e^{-y}}{ f + y + a e^{-y}}\right) $

where $\partial_t f =\dot{f}$, $\partial_x f = f'$, $d=3$ and $a$ is a parameter, now set to $a=1$. Finally, $v(d)$ is a dimension dependent parameter, to be defined later. We need two initial conditions to solve the resulting equation but the two are not independent, if we set $u'(0) = m^2$, then $u(0)= - \frac{v(d)}{d} L_{d}^{1}(m^2)$. We are not done yet, because we need to vary $m^2$ to find the true solution. For a general value of $m^2$, the solution $u^*(x)$ is going to diverge at a finite value of $x$. As $m^2$ gets closer to its 'true' value, $u^*(x)$ diverges at higher and higher $x$. I find the 'true' values of $m^2$ by dichotomy (plot the solutions and change $m^2$ accordingly) and then if i found the correct $m^2$ say, for up to 4 digits of precision, then i can look for the eigenvalues. I have two questions:

• Is there an elegant way to find the correct $m^2$?
• Is there a more efficient way to do this, than my implementation? Or maybe one that can easily be generalized for coupled PDE-s.

Clear[intL, v];
v[d_] := (2^(d + 1) \[Pi]^(d/2) Gamma[d/2])^-1
intL[d_, n_, \[Alpha]_?NumericQ][func_?NumericQ] := 
 intL[d, n, \[Alpha]][
   func] = -2 \[Alpha] NIntegrate[
    y^(d/2 - 1) ( E^-y (1 + y)  (y + E^-y \[Alpha] + func)^-n), {y, 
     0, \[Infinity]}, 
    Method -> {Automatic, "SymbolicProcessing" -> False}]
statEQ[d_, \[Alpha]_][r_] := -d u[r] + (d - 2) r u'[r] - 
  v[d] intL[d, 1, \[Alpha]][u'[r] + 2 r  u''[r]]
Clear[rStart, rEnd, bc, psol, m2, alpha];
rStart = 10^-6; rEnd = 2; alpha = 1; dim = 3;
bc = {u[rStart] == -(1/(24 \[Pi]^2)) intL[dim, 1, alpha][m2], 
   u'[rStart] == m2};
psol = ParametricNDSolveValue[{statEQ[dim, alpha][\[Rho]] == 0, bc}, 
   u, {\[Rho], rStart, rEnd}, {m2}, 
   Method -> {"EquationSimplification" -> "Residual", 
     "ParametricCaching" -> None, "ParametricSensitivity" -> None}];

(I can't directly give the initial conditions at x=0, so I chose a small rStart) Findig a good enough $m^2$ by dichotomy:

mIni = -0.2648;
fsol = psol[mIni]; // AbsoluteTiming
(*Plot[Evaluate[fsol[rho]],{rho,rStart,rEnd},PlotRange\[Rule]All]*)

Finally, obtaining the eigenvalues:

(Derivative[1]@intL[d_, n_, \[Alpha]_])[
  func_] := (Derivative[1]@intL[d, n, \[Alpha]])[
   func] = -n intL[d, n + 1, \[Alpha]][func]
linearizedDE = 
  Coefficient[
    Series[statEQ[dim, alpha][\[Rho]] /. {u[\[Rho]_] :> 
        u[\[Rho]] + \[Epsilon] \[Delta]u[\[Rho]], (Derivative[n_]@
           u)[\[Rho]_] :> (Derivative[n]@
            u)[\[Rho]] + \[Epsilon] (Derivative[
              n]@\[Delta]u)[\[Rho]]}, {\[Epsilon], 0, 
      1}], \[Epsilon]] /. u -> fsol;
linEnd = 1;
{om, nu} = 
   NDEigenvalues[
    linearizedDE, \[Delta]u[\[Rho]], {\[Rho], rStart, linEnd}, 2, 
    Method -> {"SpatialDiscretization" -> {"FiniteElement", 
        "MeshOptions" -> {"MaxCellMeasure" -> (linEnd - rStart)/
           10^3 }}}]; // AbsoluteTiming 

Answer 1

There is exact solution of stationary problem u=-v[3]/3 intL[3, 1, 1][0]=0.0102017, therefore we can consider eigensystem around this solution as follows

Clear[intL, v];
v[d_] := (2^(d + 1) \[Pi]^(d/2) Gamma[d/2])^-1
intL[d_, n_, \[Alpha]_?NumericQ][func_?NumericQ] := 
 intL[d, n, \[Alpha]][
   func] = -2 \[Alpha] NIntegrate[
    y^(d/2 - 1) (E^-y (1 + y) (y + E^-y \[Alpha] + func)^-n), {y, 
     0, \[Infinity]}, 
    Method -> {Automatic, "SymbolicProcessing" -> False}]
statEQ[d_, \[Alpha]_][r_] := -d u[r] + (d - 2) r u'[r] - 
  v[d] intL[d, 1, \[Alpha]][u'[r] + 2 r u''[r]]
Clear[rStart, rEnd, bc, psol, m2, alpha];
rStart = 10^-6; rEnd = 10; alpha = 1; dim = 3;
bc = {u[rStart] == 0.010201726677643571, u'[rStart] == 0};
sol = NDSolveValue[{statEQ[dim, alpha][\[Rho]] == 0, bc}, 
   u, {\[Rho], rStart, rEnd}];
 (Derivative[1]@intL[d_, n_, \[Alpha]_])[
  func_] := (Derivative[1]@intL[d, n, \[Alpha]])[
   func] = -n intL[d, n + 1, \[Alpha]][func]
linearizedDE = 
  Coefficient[
    Series[statEQ[dim, alpha][\[Rho]] /. {u[\[Rho]_] :> 
        u[\[Rho]] + \[Epsilon] \[Delta]u[\[Rho]], (Derivative[n_]@
           u)[\[Rho]_] :> (Derivative[n]@
            u)[\[Rho]] + \[Epsilon] (Derivative[
              n]@\[Delta]u)[\[Rho]]}, {\[Epsilon], 0, 
      1}], \[Epsilon]] /. u -> sol;
linEnd = 10;
{om, nu} = 
   NDEigensystem[{linearizedDE, 
     DirichletCondition[\[Delta]u[\[Rho]] == 0, 
      True]}, \[Delta]u[\[Rho]], {\[Rho], rStart, linEnd}, 
    10]; // AbsoluteTiming

Note, that we can take rEnd, linEnd arbitrary in this case. Visualization

Table[Plot[Evaluate[ReIm[nu[[i]]]], {\[Rho], rStart, linEnd}, 
  PlotLabel -> om[[i]], PlotRange -> All], {i, 10}]

Also we can compute solution of bouncing type. Since intL[3,1,1][f] has singular point at f=-1 we can suppose that there is a critical point solution we detect with WhenEvent[] as follows

v[d_] := (2^(d + 1) \[Pi]^(d/2) Gamma[d/2])^-1
intL[d_, n_, \[Alpha]_?NumericQ][func_?NumericQ] := 
 intL[d, n, \[Alpha]][
   func] = -2 \[Alpha] NIntegrate[
     y^(d/2 - 1) (E^-y (1 + y) (y + E^-y \[Alpha] + func)^-n), {y, 
      0, \[Infinity]}, 
     Method -> {Automatic, "SymbolicProcessing" -> False}] // Quiet
statEQ[d_, \[Alpha]_][r_] := -d u0[r] + (d - 2) r u0'[r] - 
  v[d] intL[d, 1, \[Alpha]][u0'[r] + 2 r u0''[r]]
rStart = 0; rEnd = 3.97; alpha = 1; d = 3; dim = 3;
r0 = 10^-6; m = -0.2648; eq = {-d u0[r] + (d - 2) r u0'[r] - 
    v[d] intL[d, 1, 1][u0'[r] + 2 r u0''[r]] == 0}; nds = 
 NDSolveValue[{eq, 
   u0'[r0] == m, -d u0[r0] - v[d] intL[d, 1, 1][u0'[r0]] == 0, 
   WhenEvent[u0[r] == 0, {rc = r, u0'[r] -> -u0'[r]}]}, 
  u0, {r, r0, rEnd}, 
  Method -> {"EquationSimplification" -> "Residual"}, 
  MaxSteps -> Infinity]

We plot this solution and f=u0'[r] + 2 r u0''[r] to show that solution passes point f=-1 and problem becomes stuff at r > rEnd

{Plot[nds[r],{r, r0, rEnd}, PlotRange -> All],
Plot[nds'[r] + 2 r nds''[r], {r, r0, rEnd}, PlotRange -> All]} 

Now we have position of the critical point rc=1.03115 therefore we can compute eigensystem in a range r0 <= r <=rc

(Derivative[1]@intL[d_, n_, \[Alpha]_])[
  func_] := (Derivative[1]@intL[d, n, \[Alpha]])[
   func] = -n intL[d, n + 1, \[Alpha]][func]
linearizedDE = 
  Coefficient[
    Series[statEQ[dim, alpha][\[Rho]] /. {u0[\[Rho]_] :> 
        u0[\[Rho]] + \[Epsilon] \[Delta]u[\[Rho]], (Derivative[n_]@
           u0)[\[Rho]_] :> (Derivative[n]@
            u0)[\[Rho]] + \[Epsilon] (Derivative[
              n]@\[Delta]u)[\[Rho]]}, {\[Epsilon], 0, 
      1}], \[Epsilon]] /. u0 -> nds;
linEnd = rc; NDEigenvalues[linearizedDE, \[Delta]u[\[Rho]], {\[Rho], rStart, 
  linEnd}, 2, 
 Method -> {"SpatialDiscretization" -> {"FiniteElement", 
     "MeshOptions" -> {"MaxCellMeasure" -> (linEnd - rStart)/10^2}}}]

Finally we got eigenvalues {0.636452, -1.52897}, and the question is how we can optimize initial value m. The answer depends on physics of this problem, for instance, we can make suggestion about minimum of potential in a critical point.

Answer 2

My idea is to use the domain option of the interpolation function obtained from the ndsolve. Here I set a very large value rEnd. rEnd is arbitrary, but has to be sufficiently large.

rootEq[\[Alpha]_, mass2_?NumericQ] := 
  Block[{paramSol, rhoend, rStart, rEnd, alpha, dim, fsol, delta},
   delta = 10^-6;
   rStart = delta; rEnd = 10; alpha = \[Alpha]; dim = 3;
   paramSol = 
    ParametricNDSolveValue[{statEQ[dim, alpha][\[Rho]] == 0, bc}, 
     u, {\[Rho], rStart, rEnd}, {m2},
     Method -> {"EquationSimplification" -> "Residual", 
       "ParametricCaching" -> None, "ParametricSensitivity" -> None}];
   fsol = Quiet@paramSol[mass2];
   Return[fsol["Domain"][[1, 2]]];
   ];

One then can visualize the end r*, where the solution of the differential equation becomes stiff:

dat = ParallelTable[{mass,rootEq[1,mass]}, {mass, -0.200, -0.300, -0.005}]; 
ListPlot[dat]

The x axis shows the parameter m^2, while the y axis shows the corresponding r*. As you can see, this is a quite narrow resonance curve, and $r* \to \infty$ as the corresponding $m^2_*$ is approached.

Using the function rootEq, one can simply use findroot, such as

bignumber=10;
FindRoot[rootEq[1, mass2] == bignumber, {mass2, mIni} ]

The idea is that a more precise $m^2$ corresponds to a larger $r*$. This method, with the current precision settings (basic machine precision everywhere) breaks down at around $r^*=1.8$ (corresponding to $m^2 = -0.264812$). As the workingprecison, precisiongoal and accuracygoal are increased, we can approach a higher $r^*$ using this findroot method. But now, I came across a nother problem. The eigenvalues I was looking for (and obtained from NDEigenvalues) are more sensitive to the MaxCellMeasure option, than the precision of $m^2$ beyond the 4th digit of precision... :D