Code covered by the BSD License  

Highlights from
Chebpack

image thumbnail

Chebpack

by

 

15 Jul 2011 (Updated )

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

pd.m
function [x w] = pd(n,dom,kind)
%   x: column vector of n Chebyshev points x of the kind 1 or 2 
%   or Legendre points for kind = 3
%   w: row vector of the weights w for Clenshaw-Curtis or Gauss quadrature 
%   scaled in the interval dom
% 
%   Inspired by chebpts.m from chebfun package
%   Copyright  2011, 
%   The Chancellor, Masters and Scholars of the University of Oxford, 
%   and the Chebfun Developers. All rights reserved.
%   See http://www.maths.ox.ac.uk/chebfun/ for Chebfun information.
%
%   [1] Jrg Waldvogel, "Fast construction of the Fejr and Clenshaw-Curtis
%   quadrature rules", BIT Numerical Mathematics 43 (1), pp 1--18 (2004)
%   or BIT Numerical Mathematics 46 (1), pp 195--202 (2006).
%
m = n-1;
    if kind == 1
        x = sin(pi*(-m:2:m)/(2*m+2)).';
        if nargout > 1
            w = weights1(n);
        end
    elseif kind==2 %else
        x = sin(pi*(-m:2:m)/(2*m)).';
        if nargout > 1
            w = weights2(n);
        end
    else 
        [x,w]=gauss(n);
    end
dab = diff(dom);
x = x/2*dab + (dom(1) + dab/2);
if nargout > 1, w = dab*w/2;end
function w = weights2(n) % 2nd kind
% Jrg Waldvogel, "Fast construction of the Fejr and Clenshaw-Curtis 
% quadrature rules", BIT Numerical Mathematics 46 (1), pp. 195-202 (2006).
m = n-1;  
c = zeros(1,n);
c(1:2:n) = 2./[1 1-(2:2:m).^2 ]; 
f = real(ifft([c(1:n) c(m:-1:2)]));
w = [f(1) 2*f(2:m) f(n)];
end
function w = weights1(n) % 1st kind
% Jrg Waldvogel, "Fast construction of the Fejr and Clenshaw-Curtis 
% quadrature rules", BIT Numerical Mathematics 46 (1), pp. 195-202 (2006).
l = floor(n/2)+1;
K = 0:n-l;   
v = [2*exp(1i*pi*K/n)./(1-4*K.^2)  zeros(1,l)];
w = real(ifft(v(1:n) + conj(v(n+1:-1:2))));
end
function [x,w]=gauss(n)
% Nodes and weights for Gauss quadrature on the domain dom
%
% from N. L. Trefethen, Spectral Methods in Matlab, SIAM 2000
%      http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/spectral.html, gauss.m
%       
% input:   n - number of gridpoints
%          dom - working interval
% output:  x - the (Legendre) grid (n x 1)
%          w - the weights (1 x n)
beta = .5./sqrt(1-(2*(1:n-1)).^(-2));
T = diag(beta,1) + diag(beta,-1);
[V,D] = eig(T);
x = diag(D); [x,i] = sort(x);
%xleg=x*(dom(2)-dom(1))/2+(dom(2)+dom(1))/2;
w = 2*V(1,i).^2;
end
end

Contact us