Sturm-Liouville problems

The basic Sturm-Liouville problem is

-(∂/(∂x))(p y')+q y=λ y

with boundary conditions

({{(Subscript[α, 1])y(a)+(Subscript[α, 2])y'(a), =, 0}, {(Subscript[β, 1])y(b)+(Subscript[β, 2])y'(b), =, 0}}).

If p, p', and q are continuous on [a,b] and p is never zero in [a,b], then the problem is called regular.  Otherwise, it is singular.  A regular Sturm-Liouville problem is guaranteed to have a complete, orthonormal set of eigenfunctions, while a singular problem might not.

In the most basic example, p(x)=1 and q(x)=(Subscript[α, 2])=(Subscript[β, 2])=0, so the eigenfunctions are linear combinations of trigonometric functions; this yields Fourier's theory.  The purpose of this notebook is to illustrate the existence of other examples.

Assumptions

$Assumptions={m∈Integers, m>0,n∈Integers, n>0};

Example 1

Consider the Sturm-Liouville problem

-y''=λ y.

with boundary conditions y(0)=y'(1)=0, perhaps the simplest modification of our most basic eigenvalues.

The general solution

To find the eigenvalues and eigenfunctions, the first step is to find the general solution.

Clear[y];
{genSol}=DSolve[-y''[x]==λ*y[x],y[x],x];
y[x_]=y[x]/.genSol

Finding λ

We use the boundary values to find possible values for λ.  Now, we need y'(0)=0, which implies that (Subscript[c, 1])=0.

y[0]==0

Next, we need y'(1)=0, already knowing that (Subscript[c, 1])=0.

y'[1]==0 /. C[1]->0

The only way this happens is if (Sqrt[λ])=(2n-1)π/2 for some n.  Thus our solutions to the SLP have the form

y(x)=(Subscript[c, n])sin((2n-1)π x/2)

Using the determinant

If y(x)=(Subscript[c, 1])cos((Sqrt[λ])x)+(Subscript[c, 2])sin((Sqrt[λ])x), then y'(x)=-(Subscript[c, 1])(Sqrt[λ])sin((Sqrt[λ])x)+(Subscript[c, 2])(Sqrt[λ])cos((Sqrt[λ])x) so the boundary conditions y(0)=y'(1)=0 can be written

(Subscript[c, 1])+0=0
-(Sqrt[λ])sin((Sqrt[λ]))(Subscript[c, 1])+(Subscript[c, 2])(Sqrt[λ])cos((Sqrt[λ]))=0

or, as the matrix equation

(({{1, 0}, {-(Sqrt[λ])sin((Sqrt[λ])), (Sqrt[λ])cos((Sqrt[λ]))}}))(({{(Subscript[c, 1])}, {(Subscript[c, 2])}}))=(({{0}, {0}})).

This has a non-trivial solution iff the determinant of the matrix is zero.  The determinant is just (Sqrt[λ])cos((Sqrt[λ])), which puts us where we were before.

Checking orthogonality

Clear[y];
y[n_,x_] = Sin[(2n-1)Pi*x/2];
Table[Integrate[
Sin[(2n-1)Pi*x/2]*Sin[(2m-1)Pi*x/2],{x,0,1}],
{n,1,6},{m,1,6}]//MatrixForm

Note that the 1/2s indicate that we’ll have to normalize to use in approximations.

A PDE

Let’s put this in the context of a PDE problem.  Suppose we’d like to solve (Subscript[u, t])=(Subscript[u, x x]), subject to u(0,t)=u'(1,t)=0 and u(x,0)=x(2-x).

c[n_] = 2Integrate[x(2-x)y[n,x],{x,0,1}]

Manipulate[
Plot[{x(2-x),Sum[c[n]*Exp[-t((((2n-1)Pi/2))^2)]y[n,x],{n,1,numTerms}]},
{x,0,1},PlotRange -> {0,1.2}],
{t,0,1}, {numTerms,1,10,1}]

Example II

Consider the Sturm-Liouville problem

-y''=λ y.

with boundary conditions y(0)=y(1)+y'(1)=0 - a mild variation of a text book problem.

The general solution

We certainly have the same general solution.

Clear[y];
{genSol}=DSolve[-y''[x]==λ*y[x],y[x],x];
y[x_]=y[x]/.genSol

Finding λ

We again have y(0)=(Subscript[c, 1])cos(0)+(Subscript[c, 2])sin(0)=(Subscript[c, 1])+0=0, which implies (Subscript[c, 1])=0.  
Next, we have

y(1)+y'(1)=(Subscript[c, 1])cos((Sqrt[λ]))+(Subscript[c, 2])sin((Sqrt[λ]))-(Sqrt[λ])sin((Sqrt[λ]))(Subscript[c, 1])+(Subscript[c, 2])(Sqrt[λ])cos((Sqrt[λ]))=0.

This reduces to (Subscript[c, 2])sin((Sqrt[λ]))+(Subscript[c, 2])(Sqrt[λ])cos((Sqrt[λ]))=0 or (Sqrt[λ])=-tan((Sqrt[λ])), since we already know that (Subscript[c, 1])=0.  Thus, our solutions have the specific form

y(x)=c sin((Sqrt[(Subscript[λ, n])])x)

where (Subscript[λ, n]) is one of the solutions to (Sqrt[λ])=-tan((Sqrt[λ])).  Note that the Sturm-Liousville theory guarantees the existence of infinitely many such solutions.  This equation is not solvable in closed form.  We can numerically approximate some solutions within a particular range, though.

intersections =λ/. NSolve[{-Tan[(Sqrt[λ])]==(Sqrt[λ]),0<λ<220},λ]

We can even understand the asymptotic behavior.  The following picture indicates why (Sqrt[λ]) asymptotically approaches the roots of cos(x).

tanPart=Plot[-Tan[(Sqrt[λ])],{λ,0,220}, PlotStyle -> Thickness[0.004],
Exclusions -> Cos[(Sqrt[λ])]==0, ExclusionsStyle -> Dashed,PlotRange -> {-5,15}];
sqrtPart=Plot[Sqrt[λ],{λ,0,220},
PlotStyle -> Directive[ColorData[1,2],Thickness[0.004]]];
Show[{tanPart,sqrtPart, Graphics[{PointSize[Medium],
Point[{#,(Sqrt[#])}]& /@ intersections}]}]

Checking orthogonality

Again, our solutions have the specific form

y(x)=c sin((Sqrt[(Subscript[λ, n])])x)

where (Subscript[λ, n]) is one of the solutions to (Sqrt[λ])=-tan((Sqrt[λ])) and the first few such values are

λns =λ/. NSolve[{-Tan[(Sqrt[λ])]==(Sqrt[λ]),0<λ<1000},λ]

Clear[y];
y[n_,x_] := Sin[(Sqrt[λns[[n]]])x];
Table[NIntegrate[
y[m,x]*y[n,x],{x,0,1}],
{n,1,5},{m,1,5}]//Chop//MatrixForm//Quiet

Plot[Evaluate[Table[y[n,x]/3,{n,1,5}]],{x,0,1},
AspectRatio->Automatic]

Approximating a function

Let’s use this to approximate the function f(x)=((x(1-x))^2).  Note that this function satisfies the boundary conditions y(0)=0, y(1)+y'(1)=0

Clear[n];
c[n_] :=c[n]= ((NIntegrate[x(((1-x))^2)y[n,x],{x,0,1}])/(NIntegrate[((y[n,x])^2),{x,0,1}]));

Manipulate[
Plot[{x(((1-x))^2),Sum[c[n]*y[n,x],{n,1,numTerms}]},{x,0,1},
PlotRange -> {0,0.16}],
{numTerms,1,8,1}]

Example III - Airy functions

Consider the Sturm-Liouville problem

-y''+x y=λ y

with boundary conditions y(0)=y(1)=0, perhaps the simplest Sturm-Liouville problem after Fourier's.  To find the eigenvalues and eigenfunctions, the first step is to ask Mathematica to find the general solution for us.

Remove[y];
{genSol}=DSolve[-y''[x]+x*y[x]==λ*y[x],y[x],x];
(Subscript[y, λ_])[x_]=y[x]/.genSol

The solution is expressed in terms of strange entities called Airy functions.  Of course, we should expect unfamiliar special functions, since a major point of this notebook is to encounter non-elementary Sturm-Liouville problems.

The boundary conditions state that we would like the solution to be zero at both endpoints.  Of course, we're not interested in the trivial solution, so the question is - for which values of λ are there constants, not both zero, so that the boundary conditions are simultaneously satisfied. We can write this as a matrix equation as follows.

(({{Ai(-λ), Bi(-λ)}, {Ai(1-λ), Bi(1-λ)}}))(({{(Subscript[c, 1])}, {(Subscript[c, 2])}}))=(({{0}, {0}})).

Of course, this system has a non-trivial solution precisely when the matrix

M(λ)=(({{Ai(-λ), Bi(-λ)}, {Ai(1-λ), Bi(1-λ)}}))

is singular.  Also, there is a simple test for singularity, namely the determinant must be zero.  Thus, let's define the matrix and the determinant as functions of λ and plot the value of the determinant in terms of λ.

M[λ_] = (({{AiryAi[-λ], AiryBi[-λ]}, {AiryAi[1-λ], AiryBi[1-λ]}}));
det[λ_] = Det[M[λ]]
Plot[det[λ],{λ,0,400}]

The roots are eigenvalues of the problem.  Here are all the eigenvalues less than 2000.

eigenvalues = λ/.NSolve[det[λ]==0&&0<λ<2000,λ,
WorkingPrecision->30];
N[eigenvalues]

Now, for each of the eigenvalues, there corresponds a choice of (Subscript[c, 1]) and (Subscript[c, 2]) such that (Subscript[c, 1])Ai(-λ)+(Subscript[c, 2])Bi(-λ)=0 and (Subscript[c, 1])Ai(1-λ)+(Subscript[c, 2])Bi(1-λ)=0.  Put another way, the vector ⟨(Subscript[c, 1]),(Subscript[c, 2])⟩ should map to zero under multiplication by the matrix M(λ).  Such a vector is said to lie in the null-space of the matrix. The function (Subscript[y, λ])(x)=(Subscript[c, 1])Ai(x-λ)+(Subscript[c, 2])Bi(x-λ) is then the corresponding eigenfunction.  The NullSpace command is a bit delicate numerically, which is why we requested high precision eigenvalues.  Thus, we now compute the null spaces, toss away the extra precision, define the eigenfunctions, normalize, and visualize.

nullVectors = N[First/@Table[NullSpace[M[λ]],{λ,eigenvalues}]];
eigenvalues=N[eigenvalues];
(Subscript[y, k_])[x_]:=nullVectors[[k]].{AiryAi[x-eigenvalues[[k]]],
AiryBi[x-eigenvalues[[k]]]};
Do[normSquared[k]=NIntegrate[(((Subscript[y, k])[x])^2),{x,0,1}],
{k,1,Length[eigenvalues]}]
Plot[Evaluate@Table[(Subscript[y, k])[x]/Sqrt[normSquared[k]],{k,1,8}],
{x,0,1}]

If this all worked, then the following computation should yield an identity matrix.

Table[NIntegrate[(((Subscript[y, k])[x])/(Sqrt[normSquared[k]]))(((Subscript[y, j])[x])/(Sqrt[normSquared[j]])),
{x,0,1}],{k,1,8},{j,1,8}]//Quiet//Chop//MatrixForm

OK, let's try to approximate a function.

f[x_] = (1-Abs[1-2x])/2;
Do[(Subscript[c, k])=NIntegrate[((f[x](Subscript[y, k])[x])/(normSquared[k])),{x,0,1}],
{k,1,14}];
m=Length[eigenvalues];
approx[x_] = (∑_(k=1)^m)(Subscript[c, k])(Subscript[y, k])[x];
Plot[{f[x],approx[x]},{x,0,1}]

Example IV - Bessel functions and weighted inner products

Consider the second order differential equation

-(x^2) y''-x y'=λ (x^2)y.

In this form, it is known as a Bessel equation.  While not explicitly a Sturm-Liouville equation, it is equivalent to

-(∂/(∂x))(x y')=λ x y.

This is easy to verify by expanding out the second equation using the product rule, then multiplying through by x.  (In fact, any second order, linear ODE can be put in Sturm-Liouville form using the technique of integrating factors.)  Note the x on the right hand side.  This means that the problem is not quite in the Sturm-Liouville form as presented in our text; the equation at the top of this notebook is based on the text, as well.  The slightly more general Sturm-Liouville form required here is

-(∂/(∂x))(p y')+q y=λ w(x) y.

The function w(x) is called a weight and the eigenfunctions generated by this problem will be orthogonal with respect to integration against this weight, i.e. our inner product will be defined by

(f,g)=(∫_a^b)f(x)g(x)w(x)dx.

Now, let us associate boundary conditions y'(0)=y(1)=0 with Bessel's equation.  Note the Neumann condition on the left.  Emulating the procedure from the last example, we can examine the general solution.

Remove[y];
DSolve[-(x^2)*y''[x]-x y'[x]==λ*(x^2)*y[x],y[x],x]

The solution is expressed in terms of Bessel functions.  If you try the matrix approach as before, you'll run into difficulties.  We can graph the main parts of the solution to see why.

Plot[{BesselJ[0,12x],BesselY[0,12x]},{x,0,1}]

Note that one of them (the second BesselY function, as it turns out) has an asymptote at the origin.  The other BesselJ function appears to already have derivative zero at the origin.  Thus, using these functions will boil down to finding λ such that J(0,(Sqrt[λ])1)=0.  We can do this using NSolve, as before.

eigenvalues=λ/.NSolve[BesselJ[0,Sqrt[λ]]==0 && 0<λ<1000,λ]

With the eigenvalues, we can find the eigenfunctions, normalize, and visualize.

eigenfunctions = Table[BesselJ[0,Sqrt[λ]x],{λ,eigenvalues}];
normSquares = Table[Integrate[(f^2)x,{x,0,1}],{f,eigenfunctions}];
normedEigenfunctions = 
Inner[Times,eigenfunctions,1/Sqrt[normSquares],List];
Plot[normedEigenfunctions,{x,0,1}]

Again, we should check orthonormality.

Table[NIntegrate[normedEigenfunctions[[k]]normedEigenfunctions[[j]]x,
{x,0,1}],{k,1,8},{j,1,8}]//Quiet//Chop//MatrixForm

Finally, let's try to approximate a function.

f[x_] = 1-x;
m=Length[eigenvalues];
Do[(Subscript[c, k])=NIntegrate[f[x]normedEigenfunctions[[k]]x,{x,0,1}],
{k,1,m}];
approx[x_] = (∑_(k=1)^m)(Subscript[c, k])normedEigenfunctions[[k]];
Plot[{f[x],approx[x]},{x,0,1},
AspectRatio->Automatic]

Example V - Numerical

Of course, there are many SLPs for which the general solution cannot be found in terms of any functions, special or otherwise.  Here's an example.

Remove[y];
DSolve[-D[Exp[x]y'[x],x]+Sin[x]y[x]==λ*(x^2)y[x],y[x],x]

In this case, we've got to resort to purely numerical techniques.  Let's examine this example a bit further, just to explore the possibilities.

The major difficulty is the determination of the eigenvalues, since NDSolve expects numeric quantities, not unknown parameters.  Another difficulty is the fact that we’re dealing with a boundary value problem, rather than an initial value problem.  One approach to dealing with this sort of thing is called the shooting method.  Let’s illustrate with the problem at hand.  Thus, we consider the differential equation

-(e^x) y''(x)+y(x) sin(x)==λ (x^2)y(x)

and ask, for what values of λ can the boundary conditions y(0)=y(1)=0 be satisfied.  Well, for a given λ it’s easy to solve the corresponding initial value problem with conditions y(0)=0 and y'(0)=1.  We then check the value of y(1) and adjust value of λ accordingly.  Call this process taking a shot.  Then code for a shot is as follows.

shot[λ_?NumericQ] := Module[{if,y},
if = y /. First[NDSolve[{
-D[Exp[x]y'[x],x]+Sin[x]y[x]==λ*(x^2)y[x],
y[0]==0,y'[0]==1},y,{x,0,1}]];
if[1]]

If we set λ=2, then here is the value of the solution of the IVP at x=1.

shot[2]

That doesn’t work, but shot[62] works out much better.

shot[62]

Here’s a plot of shot[λ].

pic = Plot[shot[λ],{λ,0,5000}]

As usual, the graph crosses the horizontal λ-axis a bunch of times and these roots are the eigenvalues.  Unfortunately, the NSolve command won't work on a function defined via a procedure like shot; we'll need to use FindRoot, which requires extra information regarding the approximate location of the desired roots.  Here's a fairly hacked together and probably ill-advised approach to doing so that works for this example.  The idea is to programmatically search the graph we just generated for sign changes.

pts = Cases[pic,_Line,Infinity][[1,1]];
splitPoints = Split[pts,Sign[#1[[2]]]==Sign[#2[[2]]]&];
endPoints = Flatten[{First[#],Last[#]}& /@ splitPoints,1];
endPoints = Partition[endPoints[[2;;-2]],2];
endPoints//Column

Now, each root should be bracketed by the λ-values in a pair of endpoints.  Perhaps a picture will help illustrate.

ListLinePlot[splitPoints,Epilog -> {PointSize[Large],Table[
{ColorData[1,i],Point[endPoints[[i]]]},{i,1,Length[endPoints]}]}]

Next, we use FindRoot to find the roots precisely.  Note that we're calling FindRoot in the form FindRoot[eq,{var,a,b}], which searches in the interval [a,b].  We've used the graph to verify that shot has a sign change on each such interval, so there should be a root in each one - yay, intermediate value theorem!

xEndPoints = Map[First,endPoints,{2}];
eigenvalues = λ /. Table[FindRoot[shot[λ]==0,{λ,xep[[1]],xep[[2]]}],
{xep,xEndPoints}]

Here are the corresponding eigenfunctions.

Remove[y];
k=0;
Do[k=k+1;
(Subscript[y, k])[x_] = y[x] /. First[NDSolve[{
-D[Exp[x]y'[x],x]+Sin[x]y[x]==λ*(x^2)y[x],
y[0]==0,y'[0]==1},y[x],{x,0,1}]];
(Subscript[y, k])[x_] = (Subscript[y, k])[x]/Sqrt[NIntegrate[(((Subscript[y, k])[x])^2)(x^2),{x,0,1}]],
{λ,eigenvalues}];
Plot[Evaluate[Table[(Subscript[y, k])[x],{k,Length[eigenvalues]}]],{x,0,1}]

The whole process seems complicated enough that it might be reasonable to run a check for orthogonality.  Recall the original equation:

-(e^x) y''(x)+y(x) sin(x)==λ (x^2)y(x).

The (x^2) term on the right should appear as a weight in our inner product.  Thus our ortho-check looks as follows.

Chop[Table[NIntegrate[(Subscript[y, k])[x](Subscript[y, j])[x](x^2),{x,0,1}],
{k,1,Length[eigenvalues]},{j,1,Length[eigenvalues]}],(10^(-6))]//MatrixForm