clear
echo on
%This example demonstrates the solution of the
%Sturm Liouville eigenvalue problem via conversion
%of the differential equation into a set of
%algebraic equations.  The example uses the attached
%routine slsolve.m which requires that you define
%the functions p, q, and w as well as boundary
%conditions in standard Sturm Liouville form.  You
%can extract the necessary information by typing
%help slsolve:
pause
help slsolve
pause
%Let's apply this to the simple equation:
%
%  y'' + lambda y = 0
%
%with boundary conditions:
%
%  y(0)=0,   y'(1)=0
%
pause
%The solution to this equation is just:
%
%  y = sin ((n-.5)*pi*x),   n = 1, 2, 3, ...
%
%with the eigenvalues:
%
%  lambda = (n-.5)^2 * pi^2
pause
%For this problem the functions p, q, and w are:
%  p = 1
%  q = 0
%  w = 1
%which are particularly simple.  These have been
%stored in the function files pfun.m, qfun.m, and
%wfun.m
%
%The boundary conditions are described by the array:
%
%  bc = [1,0,0,1]
pause
%Let's work the example:
bc = [1,0,0,1];
n=50;
[lambda,eigenvec]=slsolve('pfun','qfun','wfun',bc,n);
%and that's all there is to it!
pause
%Let's examine these results.  First we list the first
%five eigenvalues:
lambda(1:5).^.5/pi
pause
%Note that the eigenvalues match what we expected pretty
%well, except that the error is increasing as we go to
%higher eigenvalues.  This is clearer if we plot the 
%calculated eigenvalues against the exact values:
pause
i=[1:length(lambda)];
figure(1)
plot(i,i-.5,i,lambda.^.5/pi,'o')
%Note that the values tend to fall off very badly for
%the largest eigenvalues!
pause
%Next we look at the eigenfunctions themselves.  We shall
%only examine the first five:
x=[0:1/n:1];
figure(2)
plot(x,eigenvec(:,1:5))
%again we see that the solution starts to get a bit ragged
%for the higher eigenvectors.  Note that the leading eigen-
%vector has no nodes (zero crossings) in the interior of
%the domain, while the second has one, the third has two, and
%so forth.  This is typical of eigenvalue problems.
pause
%In conclusion, you may find this simple routine useful
%in dealing with Sturm-Liouville eigenvalue problems.  These
%differential equations are very common in transport problems
%such as the transient heating of a ball, cylinder, or slab,
%or the approach to an asymptotic velocity distribution in 
%startup flow in a pipe.  It even arises in determining the
%frequency of a violin string!
echo off


_______________________


function [lambda,eigenvec]=slsolve(p,q,w,bc,n)
%This function solves the Sturm-Liouville eigenvalue
%problem given by:
%
%  [p(x) y']' -q(x) y + lambda w(x) y = 0
%
%subject to the boundary conditions:
%
%  bc(1) y(0) + bc(2) y'(0) = 0
%  bc(3) y(1) + bc(4) y'(1) = 0
%
%over the domain 0< x <1.
%
%The function is called by the command:
%
%[lambda,eigenvec]=slsolve('pfun','qfun','wfun',bc,n);
%
%The function call requires that you provide the function
%names for the functions p, q, and w, as well as the boundary
%coefficients in the array bc.  The parameter n is the
%degree of discretization over the domain, e.g. h=1/n.
%The matrices A and W which are generated to solve the
%problem will thus be of size (n+1,n+1).  The function
%returns the eigenvalues in the array lambda (sorted by
%size in ascending order) and the matrix eigenvec which
%contains the corresponding eigenfunctions.  The eigenfunctions
%are all normalized so that the integral of the square times
%the weight function over the domain [0,1] is unity.
%
%A last note on error:  The code uses second order derivative
%approximations, so the error in the eigenvalues and eigen-
%vectors will be of O(1/n^2).  In general, the first few
%eigenvalues will be reliable, but the accuracy will deteriorate
%as you look at the higher eigenvalues, with the last few
%being meaningless.

h=1/n; %set discretization
x=[0:h:1]; %this is the array of x values

%Now we set up the arrays used in making the matrix A:
pp=zeros(1,n+1);
pm=zeros(1,n+1);
ww=zeros(1,n+1);
qq=zeros(1,n+1);
for i=2:n,
 pp(i)=feval(p,x(i)+h/2);
 pm(i)=feval(p,x(i)-h/2);
 ww(i)=feval(w,x(i));
 qq(i)=feval(q,x(i));
end

%The matrix W is easy:
weight=-diag(ww);

%The matrix A is a bit more complex.  First we do
%the main diagonal:
a=diag(-pp-pm-qq*h^2);
%and then the super and sub diagonals:
a=a+diag(pp(1:n),1);
a=a+diag(pm(2:n+1),-1);

%And now for the boundary conditions.  First at x=0:
a(1,1)=bc(1)-bc(2)*1.5/h;
a(1,2)=bc(2)*2/h;
a(1,3)=-bc(2)/2/h;

%and at x=1:
a(n+1,n+1)=bc(3)+bc(4)*1.5/h;
a(n+1,n)=-bc(4)*2/h;
a(n+1,n-1)=bc(4)/2/h;

%Finally, we divide by h^2:
a=a/h^2;

%Now we are ready to calculate the eigenvalues:
[v,d]=eig(a,weight);

%The number of eigenvalues and vectors will be less
%than the size of A and W, thus:
evals=diag(d);
i=find(isfinite(evals));
evals=evals(i);
evecs=v(:,i);

%Now we sort the eigenvalues and eigenvectors according
%to the size of the eigenvalues:
[temp,i]=sort(abs(real(evals)));
lambda=evals(i);
evecs=evecs(:,i);

%and finally, we normalize the eigenvectors so that the
%integral over the domain of the square of the eigenvector
%times the weight function equals unity. We use trapezoidal
%rule integration:
tw=diag(feval(w,x))/n;
tw(1,1)=tw(1,1)/2;
tw(n+1,n+1)=tw(n+1,n+1)/2;
enorm=abs(diag(evecs'*tw*evecs)).^.5;
for j=1:length(lambda),
 eigenvec(:,j)=evecs(:,j)/enorm(j)/sign(evecs(2,j));
end



_________________________


function y=pfun(x)
%This is the p function in the Sturm-Liouville problem:
y=ones(size(x));


_________________________


function y=qfun(x)
%This is the q function in the Sturm-Liouville problem:
y=zeros(size(x));


_________________________


function y=wfun(x)
%This is the w function in the Sturm-Liouville problem:
y=ones(size(x));