%% 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$.
%%
clear all;close 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')

