7
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Yes. It should be $y$ insted of $y^2$. I correct my typo. $\endgroup$ – Peter Dec 7 '17 at 4:23
8
$\begingroup$

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$.

$\endgroup$
  • $\begingroup$ I have Mathematica 9. It has not DEigensystem and NDEigensystem. Do these functions appear in later versions? $\endgroup$ – Peter Dec 7 '17 at 4:05
  • $\begingroup$ Can this be extended to say a fourth order system? $\endgroup$ – KraZug Dec 7 '17 at 10:21
  • $\begingroup$ @Peter, Sorry, they were introduced in V10.3. Related Q&A, with some pre-V10 solutions: (56305), (84607), (128239). $\endgroup$ – Michael E2 Dec 7 '17 at 12:08
  • $\begingroup$ @KraZug I can't get it work for y''''[x]. FEM is limited to second-order, which it will report in an error message. Otherwise, I can find no documentation about order limitations. $\endgroup$ – Michael E2 Dec 7 '17 at 12:10
  • $\begingroup$ Yes, that is my belief, it is not useful for higher order problems. $\endgroup$ – KraZug Dec 7 '17 at 13:22
7
$\begingroup$

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}]

enter image description here

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.

$\endgroup$
1
$\begingroup$

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) *)
$\endgroup$
  • $\begingroup$ It is very strange. What version of Mathematica do you use? I have Mathematica 9. After evaluating your code: {{y[x] -> 0}}. In fact, I need a numerical computation, because my differential operator has not a analytical solution. $\endgroup$ – Peter Dec 7 '17 at 4:16
  • $\begingroup$ OK, i added the method for older MMA versions. $\endgroup$ – Bill Watts Dec 7 '17 at 5:27
1
$\begingroup$

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].

$\endgroup$
0
$\begingroup$

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}]
$\endgroup$
  • $\begingroup$ I have another boundary conditions: $y[0]==1$ and $y[1]==1$. If I change the boundary condition in your programm, this code does not work. $\endgroup$ – Peter Dec 7 '17 at 4:20
  • $\begingroup$ If You change y[0] ==1 and Plot[psol[s][1]-1, {s, 0, 3.1}] val = Table[FindRoot[psol[s][1]-1, {s, s0}], {s0, {0.9, 2.9}}] $\endgroup$ – Eduard Dec 7 '17 at 7:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.