Code covered by the BSD License

# Chebpack

### Damian Trif (view profile)

15 Jul 2011 (Updated )

The MATLAB package Chebpack solves specific problems for differential or integral equations.

e=OrrSomDMS(N,R)
```function e=OrrSomDMS(N,R)
%  The script file OrrSomDMS.m computes the eigenvalues of the Orr-Sommerfeld
%  equation using NxN Chebyshev differentiation matrices from DMS
%
% From: S.C. Reddy, J.A.C. Weideman 1998
% Call: e=OrrSomDMS(128,5772); e1=OrrSomDMS(128,10000)
%

i = sqrt(-1);

[~,DM] = chebdif(N+2,2);                       % Compute second derivative
D2 = DM(2:N+1,2:N+1,2);                    % Enforce Dirichlet BCs

[x,D4] = cheb4c(N+2);                          % Compute fourth derivative
I = eye(size(D4));                        % Identity matrix

A = (D4-2*D2+I)/R-2*i*I-i*diag(1-x.^2)*(D2-I); % Set up A and B matrices
B = D2-I;

e = eig(A,B);
e=sort(e);

function [x, DM] = chebdif(N, M)

%  The function [x, DM] =  chebdif(N,M) computes the differentiation
%  matrices D1, D2, ..., DM on Chebyshev nodes.
%
%  Input:
%  N:        Size of differentiation matrix.
%  M:        Number of derivatives required (integer).
%  Note:     0 < M <= N-1.
%
%  Output:
%  DM:       DM(1:N,1:N,ell) contains ell-th derivative matrix, ell=1..M.
%
%  The code implements two strategies for enhanced
%  accuracy suggested by W. Don and S. Solomonoff in
%  SIAM J. Sci. Comp. Vol. 6, pp. 1253--1268 (1994).
%  The two strategies are (a) the use of trigonometric
%  identities to avoid the computation of differences
%  x(k)-x(j) and (b) the use of the "flipping trick"
%  which is necessary since sin t can be computed to high
%  relative precision when t is small whereas sin (pi-t) cannot.
%  Note added May 2003:  It may, in fact, be slightly better not to
%  implement the strategies (a) and (b).   Please consult the following
%  paper for details:   "Spectral Differencing with a Twist", by
%  R. Baltensperger and M.R. Trummer, to appear in SIAM J. Sci. Comp.

%  J.A.C. Weideman, S.C. Reddy 1998.  Help notes modified by
%  JACW, May 2003.

I = eye(N);                          % Identity matrix.
L = logical(I);                      % Logical identity matrix.

n1 = floor(N/2); n2  = ceil(N/2);     % Indices used for flipping trick.

k = [0:N-1]';                        % Compute theta vector.
th = k*pi/(N-1);

x = sin(pi*[N-1:-2:1-N]'/(2*(N-1))); % Compute Chebyshev points.

T = repmat(th/2,1,N);
DX = 2*sin(T'+T).*sin(T'-T);          % Trigonometric identity.
DX = [DX(1:n1,:); -flipud(fliplr(DX(1:n2,:)))];   % Flipping trick.
DX(L) = ones(N,1);                       % Put 1's on the main diagonal of DX.

C = toeplitz((-1).^k);               % C is the matrix with
C(1,:) = C(1,:)*2; C(N,:) = C(N,:)*2;     % entries c(k)/c(j)
C(:,1) = C(:,1)/2; C(:,N) = C(:,N)/2;

Z = 1./DX;                           % Z contains entries 1/(x(k)-x(j))
Z(L) = zeros(N,1);                      % with zeros on the diagonal.

D = eye(N);                          % D contains diff. matrices.

for ell = 1:M
D = ell*Z.*(C.*repmat(diag(D),1,N) - D); % Off-diagonals
D(L) = -sum(D');                            % Correct main diagonal of D
DM(:,:,ell) = D;                                   % Store current D in DM
end
function [x, D4] = cheb4c(N)

%  The function [x, D4] =  cheb4c(N) computes the fourth
%  derivative matrix on Chebyshev interior points, incorporating
%  the clamped boundary conditions u(1)=u'(1)=u(-1)=u'(-1)=0.
%
%  Input:
%  N:     N-2 = Order of differentiation matrix.
%               (The interpolant has degree N+1.)
%
%  Output:
%  x:      Interior Chebyshev points (vector of length N-2)
%  D4:     Fourth derivative matrix  (size (N-2)x(N-2))
%
%  The code implements two strategies for enhanced
%  accuracy suggested by W. Don and S. Solomonoff in
%  SIAM J. Sci. Comp. Vol. 6, pp. 1253--1268 (1994).
%  The two strategies are (a) the use of trigonometric
%  identities to avoid the computation of differences
%  x(k)-x(j) and (b) the use of the "flipping trick"
%  which is necessary since sin t can be computed to high
%  relative precision when t is small whereas sin (pi-t) cannot.

%  J.A.C. Weideman, S.C. Reddy 1998.

I = eye(N-2);                   % Identity matrix.
L = logical(I);                 % Logical identity matrix.

n1 = floor(N/2-1);               % n1, n2 are indices used
n2 = ceil(N/2-1);                % for the flipping trick.

k = [1:N-2]';                   % Compute theta vector.
th = k*pi/(N-1);

x = sin(pi*[N-3:-2:3-N]'/(2*(N-1))); % Compute interior Chebyshev points.

s = [sin(th(1:n1)); flipud(sin(th(1:n2)))];   % s = sin(theta)

alpha = s.^4;                       % Compute weight function
beta1 = -4*s.^2.*x./alpha;          % and its derivatives.
beta2 =  4*(3*x.^2-1)./alpha;
beta3 = 24*x./alpha;
beta4 = 24./alpha;
B = [beta1'; beta2'; beta3'; beta4'];

T = repmat(th/2,1,N-2);
DX = 2*sin(T'+T).*sin(T'-T);     % Trigonometric identity
DX = [DX(1:n1,:); -flipud(fliplr(DX(1:n2,:)))];   % Flipping trick.
DX(L) = ones(N-2,1);                % Put 1's on the main diagonal of DX.

ss = s.^2.*(-1).^k;              % Compute the matrix with entries
S = ss(:,ones(1,N-2));          % c(k)/c(j)
C = S./S';

Z = 1./DX;                      % Z contains entries 1/(x(k)-x(j)).
Z(L) = zeros(size(x));             % with zeros on the diagonal.

X = Z';                         % X is same as Z', but with
X(L) = [];                         % diagonal entries removed.
X = reshape(X,N-3,N-2);

Y = ones(N-3,N-2);              % Initialize Y and D vectors.
D = eye(N-2);                   % Y contains matrix of cumulative sums,
% D scaled differentiation matrices.
for ell = 1:4
Y = cumsum([B(ell,:); ell*Y(1:N-3,:).*X]); % Recursion for diagonals
D = ell*Z.*(C.*repmat(diag(D),1,N-2)-D);   % Off-diagonal
D(L) = Y(N-2,:);                              % Correct the diagonal
DM(:,:,ell) = D;                                     % Store D in DM
end

D4 = DM(:,:,4);                  % Extract fourth derivative matrix```