Solve an ODE with parameters in a boundary condition

Question

Consider the ODE:

ode = y''''[x] - 2*k^2*y''[x] + k^4*y[x] == I*k*a*((2*x - x^2 - c)*(y''[x] - k^2*y[x]) + 2*y[x]);

in which a and k are parameters and c can be regarded as an eigenvalue for which the ODE has some non-trivial solutions y[x].

It is subject to 4 boundary conditions:

bc1 = y[0] == 0;
bc2 = y'[0] == 0;
bc3 = y''[1] + k^2*y[1] + 2*y[1]/(1 - c) == 0;
bc4 = 3*I*k^2*y'[1] - k*a*(1 - c)*y'[1] - I*y'''[1] == (2*k^2 - 2*k - 1/2*k^3)*y[1]/(1 - c);

Note that bc3 and bc4 include the parameters a and/or k as well as the eigenvalue c.

I have tried to solve the ode with given parameters and plot a contour of the eigenvalue (normally a complex number) over a parameter plane:

sol[k_, a_] := ParametricNDSolveValue[{ode, bc1, bc2, bc3, bc4}, y[x], {x, 0, 1}, {k, a}]

ContourPlot[Im[c] == 0 /. c -> sol[k, a] // Evaluate, {k, 0, 2}, {a, 0, 10}]

Unfortunately, it returns nothing. Apparently, the issue should relate to c -> sol[k, a], since ParametricNDSolveValue could not obtain the eigenvalues, at least in the present form. The eignvalue c cannot be found by

ParametricNDSolveValue[{ode, bc1, bc2, bc3, bc4}, {y[x], c}, {x, 0, 1}, {k, a}]

Problem: Considering the complicated boundary conditions and the known asymptotics of the eigenvalue for $k\ll 1$,

$$c\sim 2+\mathrm{i}\frac{2}{3}(\frac{4}{5}a-1)k-\mathrm{i}\frac{k^3}{6}$$

the problem may be solved with the shooting method, in which we need an initial guess (e.g., with a small k) in order to integrate from x=0 to x=1 and then the trial solution should be iterated on $c$ to satisfy the boundary conditions at $x=1$. Can anybody give me some suggestions?

Answer 1

This question is particularly interesting to me, and I have a package that may be helpful to you here.

This particular equation is:

1. Fourth order
2. Linear
3. Inhomogeneous in the independent variable
4. Contains Robin-like boundary conditions
5. Includes the eigenvalue in the boundary conditions

I don't think that NDEigensystem can handle fourth-order systems, or the Robin-like boundary conditions. But my package is able to handle these things. It constructs the "Evans Function", an analytic function whose correspond to the eigenvalues of the original system, reducing the problem to finding roots of a smooth function of one variable (although $\mathbb{C}\to\mathbb{C}$ in this case).

First, we need to install the package:

Needs["PacletManager`"]
PacletInstall["CompoundMatrixMethod", 
    "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]

Then we first need to turn the resulting ODEs into a matrix form $\mathbf{y}'=\mathbf{A} \cdot \mathbf{y}$, using the function ToMatrixSystem:

sys[k_, a_] =  With[{a = a, k = k}, 
                ToMatrixSystem[ode, {bc1, bc2, bc3, bc4}, y, {x, 0, 1}, c]]

Here I have used With to inject the values of k and a into the system, so that e.g. sys[1,1] is well-defined.

Now we can call the function Evans to calculate the Evans function (also known as the Miss-Distance function), which utilises the method of compound matrices in the background. This function is an analytic function of the potential eigenvalue ($c$) in this case; zeroes of this function correspond to eigenvalues of the original system.

We can evaluate this function at a set of particular values of $k$, $a$ and $c$. For instance, for $c=0.8, k=1,a=1$:

 Evans[0.8, sys[1, 1]]
 (* -0.0274917 + 0.00618017 I  *)

The output of the Evans function is non-zero, so this is not an eigenvalue. You can see that the Evans function is complex for this system (due to the imaginary parts in the original equations), even for a real value of $c$.

Now the problem has turned into finding a root of this Evans function, which we can do with FindRoot:

 root = FindRoot[Evans[c, sys[1, 1]], {c, 2}]
 {c -> 1.30138 - 0.0344041 I}

If we substitute this value into the Evans function, we see it is zero to machine precision:

 Evans[c /. root, sys[1, 1]]
 (* 1.99493*10^-17 + 2.49366*10^-18 I *)

This method appears to give solutions that are close to your asymptotic result for $k<<1$:

asym = (2 + (2 I/3) ((4 a/5) - 1) k - (I k^3)/6);
asym /. a -> 1 /. k -> 0.001
(* 2. - 0.000133334 I *)

c /. FindRoot[Evans[c, sys[0.001, 1]], {c, 2}]
(* 2. - 0.000132668 I *)

Note there may be multiple roots, and FindRoot will only find one at a time. It can also fail to find a root at all if you don't start "close enough" (for some definition of close).

It is often helpful to loop through a set of points, taking the previous eigenvalue as the initial guess. I'll also turn off the error message for k=0, where the eigenvalue is only present in a boundary condition, not the ODE (it is fine, I just need to edit my code to not complain about that).

Off[Evans::noEigenvalue]
cstart = 2;
pts = Monitor[
  Table[root = c /. FindRoot[Evans[c, sys[k, 1]], {c, cstart}]; 
   val = Abs[Evans[root, sys[k, 1]]]; 
   If[val > 10^-10, 
    Print["Not converged at k = ", k, " , Abs(Evans) = ", val]]; 
   cstart = root; {k, root},
   {k, 0, 5, 0.1}
   ], {c, k}]

Another way is to keep a track of all the previous successful eigenvalues, and take the Nearest of those. We start with a few values (just found quickly):

goodPoints = {{0, 1} -> 2, {0, 100} -> 2};

and we make a function that will start at the nearest best guess. If the magnitude of the variation in the two axes are very different this won't work very well without modification, but it seems to be fine here. Note that we check whether FindRoot has converged sufficiently, if not then we don't output it as a point. You can also add additional steps to take in that case (such as trying some other starting values for instance).

Clear[fr]
fr[k_?NumericQ, a_?NumericQ] := Module[{root, val},
  root = c /. 
    FindRoot[
     Evans[c, sys[k, a]], {c, First@Nearest[goodPoints, {k, a}]}];
  val = Abs[Evans[root, sys[k, a]]];
  If[val < 10^-8, AppendTo[goodPoints, {k, a} -> root]; fr[k,a] ={k, a, root}, 
   Sequence[]]
  ]

pts = Monitor[Table[fr[k, a], {k, 0, 3, 0.1}, {a, 0, 10, 1}], {k, a, c}];

And then you can use ListPlot3D (or ListPointPlot3D) to plot the 3D surface of how $c$ changes with $a$ and $k$:

ListPlot3D[{#[[1]], #[[2]], Re[#[[3]]]} & /@ Flatten[pts, 1], 
  PlotRange -> All, AxesLabel -> {"k", "a", "Re(c)"}, 
  LabelStyle -> (FontSize -> 16)]
ListPlot3D[{#[[1]], #[[2]], Im[#[[3]]]} & /@ Flatten[pts, 1], 
  PlotRange -> All, AxesLabel -> {"k", "a", "Im(c)"}, 
  LabelStyle -> (FontSize -> 16), ImageSize -> 600]

For the eigenfunctions, I don't have any code currently in my package to get those out. But the method that I usually use for non-stiff equations is to find the eigenvalue as above, and then integrate the original equations using NDSolve, replacing one boundary condition with an arbitrary value. For example, for $k=1,a=1$, we first find the root:

root = FindRoot[Evans[c, sys[1, 1]], {c, 2}];

Then we integrate, but instead of using the fourth boundary condition I replace it with $y''(0)=1$. This works because eigenfunctions are unique only up to a constant, so we can set this arbitrarily.

sol = NDSolve[{ode, bc1, bc2, bc3, y''[0] == 1} /. root /. k -> 1 /. 
    a -> 1, y, {x, 0, 1}][[1]]

Plotting the real and imaginary parts:

Plot[Evaluate@ReIm[y[x] /. sol], {x, 0, 1}]

We can check that the resulting solution meets the boundary conditions (to within the accuracy of the NDSolve and FindRoot, you can up those if you want to meet it more precisely):

({bc1, bc2, bc3, bc4} /. Equal -> Subtract) /. root /. k -> 1 /. 
  a -> 1 /. sol // Chop
(* {0, 0, -3.44704*10^-9 - 2.21806*10^-8 I, 2.3183*10^-8 + 1.51406*10^-7 I}

For further details on the method, as well as more examples please see my other answers on this site, or my github page. The notebook I have on github has a number of examples, including on the continuation.