Spectral problem for differential operator

Question

I want to find numerically eigenvalue and eigenfunctions some nontrivial differential operator. But I can not find, how to do it in Mathematica.

For the sake of simplicity let us discuss the trivial problem:

$ \frac{d^2 y}{d x^2}=-\omega^2 y \quad y(0)=y(1)=0$

How to find solutions of this trivial problem using Mathematica?

Answer 1

In versions 10.3 and higher, there's DEigensystem and NDEigensystem:

DEigensystem[{-Laplacian[y[x], {x}], 
  DirichletCondition[y[x] == 0, True]}, y[x], {x, 0, 1}, 4]
(*
  {{π^2, 4 π^2, 9 π^2, 16 π^2},                      <-- eigenvalues
   {Sin[π x], Sin[2 π x], Sin[3 π x], Sin[4 π x]}}   <-- eigenfunctions
*)

They require a linear operator, but it doesn't have to be trivial. The above is for $-{d^2y\over dx^2} = \lambda\,y$.

Answer 2

The answers you have so far work nicely for the trivial problem that you posed, but fail when you try a more complex problem, where you need to do everything numerically.

I asked a very similar question to this for solving a 10th order BVP eigenvalue problem. I wrote some code in order to solve it using the Compound Matrix Method to calculate the Evans function - an analytic function of $\omega$ whose roots correspond to eigenvalues of the original problem, which I've just put on GitHub.

Install the package via:

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

Then for your example, we first need to convert the system of ODEs and BCs into the correct matrix form,

Needs["CompoundMatrixMethod`"]
sys = ToMatrixSystem[y''[x] == - ω^2 y[x], {y[0] == 0, y[1] == 0}, 
      y, {x, 0, 1}, ω]

Now Evans[ω, sys] gives a function that can be evaluated at any complex $\omega$. Roots of this function give the eigenvalues of the original BVP.

Plot[Evans[ω, sys], {ω, 0, 15}]

Here you can see the roots are given by $\omega = n \pi$. The method works for a general system of equations, and can do infinite or semi-infinite domains.

Answer 3

MMa can't do the y^2, but if you mean this:

    sol = DSolve[{y''[x] + \[Omega]^2 y[x] == 0, y[0] == 0, y[1] == 0},y[x], x]

(* y[x]->C[1] Sin[x Sqrt[\[Omega]^2]]n \[Element] Integers && n >= 1 && \[Omega]^2==n^2 \[Pi]^2 *)

If you really want to do it with the y^2 and want an analytical solution, MMa can solve the ode without the conditions. Have to use FindRoot to solve for the constants, since the equations are trancendental.

The above eigenvalues and eigenvectors use M11. For older versions do the following. Similar to what you would do on paper.

pde = y''[x] + \[Omega]^2 y[x] == 0

DSolve without conditions

sol = DSolve[pde, y[x], x] // Flatten
(* c1 Cos[x \[Omega]] + c2 Sin[x \[Omega]] *)

y[x_] = y[x] /. % /. {C[1] -> c1, C[2] -> c2}

bc1 y[0]==0

c1 = c1 /. Solve[y[0] == 0, c1][[1]]

(0)

Now if we just solve for c2 from the second bc we will get the trivial answer 0. Instead

y[1]==0

c2 Sin[\[Omega]] == 0

Reduce[{%, \[Omega] > 0}, \[Omega]]

(c2 == 0 && \[Omega] > 0) || (C[1] \[Element] 
    Integers && ((C[1] >= 1 && \[Omega] == 2 \[Pi] C[1]) || (C[1] >= 
        0 && \[Omega] == 2 \[Pi] C[1] + \[Pi])))

Boils down to

\[Omega] = n Pi

$Assumptions = n \[Element] Integers

y[x]
(* c2 Sin[Pi n x] *)

Can normalize if you want

c2 = c2 /. Solve[Integrate[y[x]^2, {x, 0, 1}] == 1, c2][[2]]
(* Sqrt[2] *)

y[x]
(* Sqrt[2] sin(Pi n x) *)

Answer 4

Here's a method to fudge NDSolve for the solution of eigenvalue problems: you simply pretend that \[Omega] is also a function of x and add it as an extra function to solve for. You can do this by specifying that \[Omega]''[x] == 0, and giving it boundary conditions that allow it to be any constant value.

NDSolve[{
  y''[x] == -\[Omega][x] y[x],
  \[Omega]''[x] == 0,
  y[0] == 1,
  y[1] == 1,
  \[Omega]'[0] == 0,
  \[Omega]'[1] == 0
  },
 {y, \[Omega]},
 {x, 0, 1}
 ]

Now you don't have any control over which value of \[Omega] it's going to find, but I think you can use the "Shooting" method with appropriate initial guesses to find other values of \[Omega].

Answer 5

Are You sure about the y^2?

According to documentation, it might be solved

psol = ParametricNDSolveValue[{y''[x] == -\[Omega]^2 y[x]^2, 
y[0] == 0, y'[0] == \[Omega]^2}, y, {x, 0, 1}, {\[Omega]}];
Plot[psol[s][1], {s, 0, 3.1}]
val = Table[FindRoot[psol[s][1], {s, s0}], {s0, {0.9, 2.9}}]
Plot[Evaluate[psol[s][t] /. val], {t, 0, 1}]