No BSD License  

Highlights from
ChebyshevTools

from ChebyshevTools by Michael Watson
Tool box for solving ODE/PDEs using spectral Chebyshev differentiation matrices.

g=cheb2galerkin(f,typ)
function g=cheb2galerkin(f,typ)

%this function transforms from the chebyshev basis to the galerkin basis

%   1) Dirichlet/Dirichlet Mason Galerkin
%   2) Dirichlet/Dirichlet Julien Galerkin
%   3) Dirichlet/Dirichlet Trefethen Galerkin
%   4) u=u'=0
%   5) u=u''=0
%   6) Neumann/Neumann
%   7) u'[-1]-u[-1]/-1=u'[1]-u[1]==0  %stress free angular velocity
%   8) u[1]=u[-1]=u''[1]+u'[1]-u[1]=u''[-1]-u'[-1]-u[-1]=0 % stress free velocity potential
%   9) u'[1]=u'[-1]=u'''[1]=u'''[-1]==0 %Neumann/Neumann and third deriv
%  10) u[1]=u[-1]=u''''[1]=u''''[-1]==0 %Dirichlet/Dirichlet and fourth deriv


switch typ
    case 1 %Mason
       g=f(3:end);
    case 2 %Julien
       % Solve via bi-diagonal backsolve instead of matrix inversion
       g=zeros(length(f)-2,1);
       %The bottom 2 entries are the same
       g(end)=f(end);
       g(end-1)=f(end-1);
       for j=length(g)-2:-1:1
           g(j)=f(j+2)+g(j+2);
       end
    case 3 %Trefethen
       % Solve via bi-diagonal backsolve instead of matrix inversion
       g=zeros(length(f)-2,1);
       %The bottom 2 entries are the same
       g(end)=f(end);
       g(end-1)=f(end-1);
       % directly calculate the next 2
       g(end-2)=f(end-2)+2*g(end);
       g(end-3)=f(end-3)+2*g(end-1);
       for j=length(g)-4:-1:1
           g(j)=f(j+2)+2*g(j+2)-g(j+4);
       end
    case 4 %BI harmonic Use matrix to start
        L=length(f)-4;
        D=ones(L+4,1);
        L2=[];
           for k=0:L+1;
              L2=[L2,-2*(k+2)/(k+3)]; 
           end
        L4=[];
           for k=0:L-1;
              L4=[L4,(k+1)/(k+3)]; 
           end 
        A=sparse(diag(D))+sparse(diag(L2,-2))+sparse(diag(L4,-4));
        g=A\f;
        g=g(1:end-4);
    case 5 %u[1]=u[-1]=u''[1]=u''[-1]=0
        L=length(f)-4;
        D=ones(L+4,1);
        L2=[];
           for k=0:L+1;
              L2=[L2,-(1+((k+1)*(2*k^2+4*k+3)/((k+3)*(2*k^2+12*k+19))))]; 
           end
        L4=[];
           for k=0:L-1;
              L4=[L4,((k+1)*(2*k^2+4*k+3)/((k+3)*(2*k^2+12*k+19)))]; 
           end 
        A=sparse(diag(D))+sparse(diag(L2,-2))+sparse(diag(L4,-4));
        g=A\f;
        g=g(1:end-4);
     case 6 %u'[1]=u'[-1]==0 %Neumann/Neumann
        L=length(f)-2;
        D=[1];
           for k=1:L+1;
              D=[D,(k+2)^2/(2*k^2+4*k+4)]; 
           end
        L2=[0];
           for k=1:L-1;
              L2=[L2,-((k)^2)/(2*k^2+4*k+4)]; 
           end 
        A=sparse(diag(D))+sparse(diag(L2,-2));
        g=A\f;
        g=g(1:end-2);
     case 7 %u'[-1]-u[-1]/-1=u'[1]-u[1]==0
        L=length(f)-2;
        D=[];
           for k=0:L+1;
              D=[D,((k+2)^2-1)/(2*(k+1)^2)]; 
           end
        L2=[];
           for k=0:L-1;
              L2=[L2,-((k^2-1))/(2*(k+1)^2)]; 
           end 
        A=sparse(diag(D))+sparse(diag(L2,-2));
        g=A\f;
        g=g(1:end-2);
     case 8 %u[1]=u[-1]=u''[1]+u'[1]-u[1]=u''[-1]-u'[-1]-u[-1]=0
        L=length(f)-4;
        D=ones(L+4,1);
        L2=[];
           for k=0:L+1;
              L2=[L2,-2*(k+2)*(k^2+4*k+9)/((k+3)*(k^2+6*k+11))]; 
           end
        L4=[];
           for k=0:L-1;
              L4=[L4,(k+1)*(k^2+2*k+3)/((k+3)*(k^2+6*k+11))]; 
           end 
        A=sparse(diag(D))+sparse(diag(L2,-2))+sparse(diag(L4,-4));
        g=A\f;
        g=g(1:end-4); 
     case 9 %u'[1]=u'[-1]=u'''[1]=u'''[-1]==0 %Neumann/Neumann and third deriv
        L=length(f)-4;
        D=ones(L+4,1);
        L2=[];
           for j=0:L+1;
              L2=[L2,(-2*j^2*(11 + 2*j*(4 + j)))/((2 + j)*(3 + j)*...
                 (15 + 2*j*(6 + j)))]; 
           end 
        L4=[];
           for j=0:L-1;
              L4=[L4,(j^2*(1 + j)*(-1 + 2*j*(2 + j)))/ ...
                  ((3 + j)*(4 + j)^2*(15 + 2*j*(6 + j)))]; 
           end 
        A=sparse(diag(D))+sparse(diag(L2,-2))+sparse(diag(L4,-4));
        
        g=A\f;
        g=g(1:end-4);
     case 10 %u[1]=u[-1]=u''''[1]=u''''[-1]==0 %Dirichlet/Dirichlet and fourth deriv
        L=length(f)-4;
        D=ones(L+4,1);
        L2=[];
           for j=0:L+1;
              L2=[L2,(-2*(210 + j*(4 + j)*(53 + 2*j*(4 + j))))/ ...
                     ((4 + j)*(5 + j)*(21 + 2*j*(6 + j)))]; 
           end 
        L4=[];
           for j=0:L-1;
              L4=[L4,(j*(-5 + j*(1 + 2*j*(1 + j))))/...
                     ((4 + j)*(5 + j)*(21 + 2*j*(6 + j)))]; 
           end 
        A=sparse(diag(D))+sparse(diag(L2,-2))+sparse(diag(L4,-4));
        
        g=A\f;
        g=g(1:end-4);
end

Contact us at files@mathworks.com