Nonlinear ODE eigenvalue problem

Question

How does one find eigenvalues $\lambda$ of the following problem?

$$ \frac{\mathrm{d}^2 u}{\mathrm{d}x^2} = \lambda \left( -u + u^2 \right),$$

$$ u(0) = u(1) = 0. $$

Can this be tackled by Mathematica's built-in method DEigensystem? It seems like the right-hand side of the equation must have the form $\lambda u$ only to make DEigensystem work.

Answer 1

We want to find the value of $\lambda$ for which there exists a solution of the differential equation $u'' = \lambda ( - u + u^2)$, subject to the boundary conditions $u(0) = u(1) = 0$. We can do this by a shooting algorithm, i.e., solving the related initial value problem $$ u'' = \lambda (- u + u^2), \qquad u(0) = 0, \; u'(0) = v. $$ The value of $u(1) \equiv b$ can then be seen to be a function of $\lambda$ and $v$. If we can find the value(s) of $\lambda$ and $v$ for which $b = 0$, then we have solved this "eigenvalue" problem.

To do this, we can use ParametricNDSolve, viewing $\lambda$ and $v$ as "parameters". We can then use ContourPlot to visualize the values of $\lambda$ and $v$ for which $u(1) = 0$. Here, the function boundaryval returns $u(1)$. A wrinkle occurs because this ODE, being non-linear, does not always have a solution over the full range $[0,1]$ for arbitrary initial conditions at 0; and sometimes the range of the InterpolatingFunction returned by ParametricNDSolve does not include $x = 1$. To account for this, we use Check; if the InterpolatingFunction returns the error dmval, we know that there is not a solution of this ODE on the domain $[0,1]$ for these parameters, and we return Null instead of a numerical value.

eqns = {u''[x] == λ (-u[x] + u[x]^2), u[0] == 0, u'[0] == v}
sol = ParametricNDSolve[eqns, u, {x, 0, 1}, {λ, v}]
boundaryval[λ_, v_] := Check[u[λ, v][1] /. sol, Null, {InterpolatingFunction::dmval}]
ContourPlot[boundaryval[λ, v] == 0, {λ, -100, 100}, {v, -20, 20},
  FrameLabel -> Automatic, PlotPoints -> 100]

(This plot took about 5 minutes to generate on my laptop.)

The line at $v = 0$ comes from the fact that when $v = 0$, the solution to the ODE is $u(x) = 0$ for all $x$, regardless of the value of $\lambda$. This central line therefore corresponds to the trivial solution to the equation. However, the curved lines correspond to non-zero values of $\lambda$ and $v$ that yield $u(1) = 0$; and it appears that there exists at least one non-zero solution to the original "eigenvalue" problem for any value of $\lambda$, including negative values. Note that for small values of $v$, the curves intersect the $\lambda$ axis at $\pi^2 \approx 9.9$, $4 \pi^2 \approx 39.5$, and $9\pi^2 \approx 88.8$; these correspond to solutions that remain small ($u \ll 1$) over the whole domain, and so the solutions are "nearly linear". In this limit, we recover the eigenvalues $n^2 \pi^2$ of the linear problem.

(I'm not sure whether the "dotted" lines are real solutions or spurious. I believe they correspond to the borders of regions in the $\lambda v$-plane where the initial value problem cannot be solved over the whole domain $[0,1]$.)

To see why we have discrete eigenvalues in the linear case but continuous eigenvalues in the non-linear case, we can also apply the same technique to the linear version of this equation:

lineareqns = {u''[x] == λ (-u[x]), u[0] == 0, u'[0] == v}
sol = ParametricNDSolve[lineareqns, u, {x, 0, 1}, {λ, v}]
boundaryval[λ_, v_] := Check[u[λ, v][1] /. sol, Null,{InterpolatingFunction::dmval}]
ContourPlot[boundaryval[λ, v] == 0, {λ, -100, 100}, {v, -20, 20}, FrameLabel -> Automatic, PlotPoints -> 100]

(This plot is generated much more quickly by Mathematica—about 10 seconds on my machine.)

We can see that other than the line $v = 0$, the contours for which $b(\lambda,v) = 0$ are vertical lines, and thus only discrete values of $\lambda$ are eigenvalues of this equation.

What's going on? In the linear case, suppose that $\{\lambda,v\}$ is a solution of $b(\lambda, v) = 0$. This means that there exists a solution $u(x)$ such that $u(0) = 0$, $u'(0) = v$, and $u(1) = 0$. But since the equation is linear for a fixed value of $\lambda$, $\tilde{u}(x) = \alpha u(x)$ will be a solution with $\tilde{u}(0) = 0$, $\tilde{u}'(0) = \alpha v$, and $\tilde{u}(1) = 0$, and with the same value of $\lambda$. Since this is true for any value of $\alpha$, it must be the case that varying $v$ does not change $b$ when $b = 0$. In other words, for the linear equation, $$ b(\lambda, v) = 0 \quad \Rightarrow \quad \left( \frac{\partial b}{\partial v} \right)_\lambda = 0. $$ Thus, the "curves" of constant $b$ in the $\lambda v$-plane point directly in the $v$-direction. No such argument can be made for the non-linear equation, since for a fixed value of $\lambda$, a multiple of a solution $u(x)$ to the ODE is generally not itself a solution.