Sturm-Liouville eigenvalue problem solver including projections

This script solves the classical Sturm-Liouville problem. This problem appears throughout applied mathematics and is similar to the traditional eigenvalue problem from linear algebra. However rather than a matrix we now have an operator and the variables are functions. We then project a function onto the basis formed by the eigenfunctions. The equations are:

$\frac{d}{dx} (p(x) \frac{du}{dx}) + q(x) u(x) = -\lambda \sigma(x) u(x)$

$p(x) u''(x) + p'(x) u'(x) + q(x) u(x) = -\lambda \sigma(x) u(x)$

The inner product is $<u,v>=\int_{-1}^1 u(x) v(x) \sigma(x) dx$ .

Calculate Modes
Figure 1
Figure 2
Figure 3
Projections
Figure 4

clear all
format long, format compact

Number of points and the derivative and integration material.

N=50;
[D,x]=cheb(N); D2=D^2; % Differentiation matrices.
[xi wt]=clencurt(N); % Integration coefficients.

Note the computational domain is [-1,1]

Define the various functions.

numdx=1e-8;
myp=@(x) 1+exp(-(x/0.2).^2);
mypx=@(x) (myp(x+numdx)-myp(x-numdx))/(2*numdx);
myq=@(x) x;
mysig=@(x) 1+0.5*sin(x*pi);

Build the operators.

myoplh=diag(myp(x))*D2+diag(mypx(x))*D+diag(myq(x));
myoprh=diag(mysig(x));

The boundary conditions are u=0 so just strip off the first and last row and column.

myoprh=myoprh(2:end-1,2:end-1);
myoplh=myoplh(2:end-1,2:end-1);

Solve the eigenvalue problem and sort the eigenvalues.

[myphis0 mylams0]=eig(myoplh,myoprh);
[mylams mylamsi]=sort(diag(mylams0),'descend');

Set the number of modes you wish to consider.

num_modes=20;
modesu=zeros(N+1,num_modes);
modesu(2:N,:)=myphis0(:,mylamsi(1:num_modes));
modese=zeros(size(modesu));

Calculate Modes

for k=1:num_modes

    unow=modesu(:,k);

Normalize using the inner product.

    normalizationu(k)=sqrt(sum(wt'.*unow.^2.*mysig(x)));
    modesu(:,k)=unow/normalizationu(k);

Notice the way the Clenshaw Curtis weights come in.

end

Figure 1

figure(1)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',12,...
    'DefaultTextFontWeight','bold','DefaultAxesFontSize',12,...
    'DefaultAxesFontWeight','bold');
plot(mylams(1:20),'bo')
xlabel('mode number')
ylabel('frequency')
grid on

Figure 2

figure(2)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',12,...
    'DefaultTextFontWeight','bold','DefaultAxesFontSize',12,...
    'DefaultAxesFontWeight','bold');
plot(x,myp(x),'k-',x,myq(x),'b-',x,mysig(x),'r')
xlabel('x')
xlabel('Sturm Liouville functions')
legend('p(x)','q(x)','\sigma(x)')

Figure 3

figure(3)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',12,...
    'DefaultTextFontWeight','bold','DefaultAxesFontSize',12,...
    'DefaultAxesFontWeight','bold');
for modei=1:6
    subplot(6,1,modei)
    plot(x,modesu(:,modei),'k-')
    xlabel('x')
    ylabel('e-funcs')
    grid on
    axis([-1 1 -2 2])
end

Projections

Now for some projections. First choose the function to project. Anything that satisfies the BCs is OK.

f=sin(2*pi*x);
myapp=zeros(size(f));
for k=1:num_modes/2
    phinow=modesu(:,k);
    coeff(k)=sum(wt'.*phinow.*f.*mysig(x));
    myapp=myapp+coeff(k)*phinow;
end
myapphalf=myapp;
for k=(num_modes/2+1):num_modes
    phinow=modesu(:,k);
    coeff(k)=sum(wt'.*phinow.*f.*mysig(x));
    myapp=myapp+coeff(k)*phinow;
end

Figure 4

figure(4)
clf
set(gcf,'DefaultLineLineWidth',2,'DefaultTextFontSize',12,...
    'DefaultTextFontWeight','bold','DefaultAxesFontSize',12,...
    'DefaultAxesFontWeight','bold');
subplot(2,1,1)
plot(x,f,'k',x,myapp,'b',x,myapphalf,'r')
xlabel('x')
ylabel('f and its approximation')
grid on
title(['f and its approximation using ' int2str(num_modes) ' and ' int2str(num_modes/2) ' us'])
subplot(2,1,2)
plot(x,f-myapp,'b',x,f-myapphalf,'r')
legend('full approx','half approx')
grid on
xlabel('x')
ylabel('error')