# Stationary solution of PDE and eigenvalues

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

• This is not PDE, this is integrodifferential equation. Your problem is not state clear. What is domain for solution of you equation? Is it {t,0,Infinity}, {x,0,2}? What do you try to minimize with respect to m? Mar 5, 2021 at 13:22
• In principle the domain is {t,0,Infinity}, {x,0,Infinity}. By fine tuning the parameter m^2 we can increase the range {x,0, something(m^2) }, where you can obtain a solution where the function u does not diverge. If I set a high rEnd, the ndsolve will signal, that the system became stiff. In that sense, you do not minimize anything w.r.t m^2, you need to fine tune it. Mar 5, 2021 at 13:40
• What do you expecting at x->Infinity? Is it zero or some limited value? Mar 5, 2021 at 13:45
• But intL[3, 1, 1][0] exists, therefore there is limited value at infinity u=-v[3]/3 intL[3, 1, 1][0]=0.0102017. Then we can normalize solution so that u->0 at infinity and use 2 conditions u[0]=-0.0102017 and u'[rEnd]=0 with rEend=10 for example in numerical computation. Mar 5, 2021 at 14:15
• Not denominator, but parameter func becomes close to -1. Mar 7, 2021 at 10:13

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.

• It does not seem to be correct to give u'[rStart]=0 as initial condition. The function sol in your answer is a consant function (which is indeed a stationary solution but I'm looking for the non-trivial one) and therefore the eigenvalues are not correct either. The resulting u function has to be something like fsol in my question and the first two eigenvalues om and nu are also the correct (or very close to the correct) ones in my question. Mar 5, 2021 at 18:33
• @dzsoga What do you mean under nontrivial? If solution is unique, then there is only one solution u=-v[3]/3 intL[3, 1, 1][0]=0.0102017 not diverges at x->Infinity. If solution not unique then there are many solutions non restricted at x->Infinity and therefore your problem not defined well. Hence you need to define the problem so that the solution is unique. Your own solution diverges at x->Infinity. Mar 5, 2021 at 18:48
• The function u here is a thermodynamic potential. The t variable the reduced temperature. You are right, the PDE for it has multiple stationary solutions. In fact it has two. The solutions are distinguished by their eigenvalues(critical exponents). The nontrivial u solution corresponds to a second order phase transition. Please observe, how my fsol changes as mIni is changed. Mar 5, 2021 at 20:04
• @dzsoga How did you derive this equation? Mar 5, 2021 at 20:31
• It is the Wetterich equation for a Z_2 symmetric scalar model (which happens to be in the Ising universality class) at the Local Potential Approximation of the effective action. Mar 5, 2021 at 20:39

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

• Ok! It is a nice approach (+1). Probably we can also use physical reasons to solve this problem. Mar 9, 2021 at 15:21