pd_old.m

MDPSAS工具箱是马里兰大学开发的

function [z,w,dz,ddz,Q] = pd(geom,ndisc,axlim)

% PD.M Computes quadrature weights, discretized coordinates, Lagrange-
%      interpolation polynomial-based differentiation arrays, and
%      discretzed Jacobi polynomials for a 1-dimensional domain. Call
%      as
%
%      [z,w,dz,ddz,Q] = pd(geom,ndisc,axlim);
%
%      INPUT Variables: 
%          geom  : geometry - 'slab', 'cyln', 'sphe', 'peri'
%          ndisc : number of discretization points, including endpoints
%          axlim : [lower upper] bounds on discretized interval
%
%      OUTPUT Variables:
%          z     : (ndisc by 1) set of discretization points in
%                  the specified interval;
%          w     : (ndisc by 1) quadrature weight array; note that
%                  the quadrature array contains the factor z.^a;
%          dz    : (ndisc by ndisc) first order derivative array;
%          ddz   : (ndisc by ndisc) Laplacian operator array;
%          Q     : discretized Jacobi polynomials in column format.
%
%      EXAMPLE: Generate the location array, quadrature weights,
%               differentiation arrays, and Jacobi polynomials
%               for a cylindrical physical domain. Plot the first
%               5 polynomials and demonstrate their orthogonality.
%
%               [r,w,dr,ddr,Q] = pd('disk',20);
%               plot(r,Q(:,1:4))
%               wip( Q(:,1:5),Q(:,1:5),w.*r.*(1-r) )

% Written by Raymond A. Adomaitis, Hsiao-Yung Chang, and Yi-hung Lin 
% as part of MWRtools, version 2.

% Copyright (c) by Raymond A. Adomaitis, 1998-2002

if nargin < 3, axlim = [0 1]; end

switch lower(geom)
   case 'slab', a = 0;
   case 'cyln', a = 1;
   case 'sphe', a = 2;
   case 'peri', % not relevant
   otherwise, disp('unknown geometry or not yet supported')
   end

if strcmp(lower(geom),'peri')

% For periodic domains, axlim are set to [0, 2*pi - delT]

   if ~mod(ndisc,2)
      error('Odd number of points required w/ periodic coords')
      end
   z = 2*pi*[0:ndisc-1]'/ndisc;
   w = ones(size(z))/ndisc;

   Q   = dFs(z,ndisc);
   dQ  = zeros(size(z));
   ddQ = zeros(size(z));

   for i = 1:(ndisc-1)/2
      dQ(:,2*i)    = -i*sin(i*z);
      dQ(:,2*i+1)  =  i*cos(i*z);
      ddQ(:,2*i)   = -i^2*cos(i*z);
      ddQ(:,2*i+1) = -i^2*sin(i*z);
      end

   dz  = dQ*inv(Q);
   ddz = ddQ*inv(Q);

   return
   end

% Non-periodic domain

z        = qpts(a,ndisc); 
w        = qweights(z,a,ndisc); 
[dz ddz] = lDM(z,ndisc);
Q        = dJs(z,a+1,ndisc-1,2);

% Adjust for axis limits

zBar = z;

z0  = axlim(1);
z1  = axlim(2);
z   = z0 + (z1-z0)*zBar;

dz  = dz/(z1-z0);
ddz = ddz/(z1-z0)^2;
w   = (z1-z0)*w;

% Laplacian operator for a > 0 requires asymptotic value when z0 = 0

if a ~= 0 & z0 > 0
   ddz = a*dz./z(:,ones(ndisc,1)) + ddz;
   end

if a ~= 0 & z0 == 0

   v = [2:ndisc];

   ddz(v,:) =  a*dz(v,:)./z(v,ones(ndisc,1)) + ddz(v,:);
   ddz(1,:) = (a+1)*ddz(1,:);
   end

% The quadrature weights, scaled for non-unit intervals

switch a

case 0
   w = w;
case 1
   w0 = qweights(zBar,0,ndisc); 
   w = z0*w0 + w;
   w = 2*pi*(z1-z0)*w;
case 2
   w0 = qweights(zBar,0,ndisc); 
   w1 = qweights(zBar,1,ndisc); 
   w = z0^2*w0 + 2*z0*(z1-z0)*w1 + (z1-z0)*w;
   w = 4*pi*(z1-z0)*w;

   end

% ------------- quadrature points ------------------ %

function [x] = qpts(a,N); 
  
x  = [0; rJ(a+1,N-2); 1];   

% ------------- quadrature weights ------------------ %

function [w] = qweights(x,a,N); 
  
dx = x(:,ones(1,N))' - x(:,ones(1,N)); % difference between all points
p  = [ones(1,N); cumprod(dx)];
dp = zeros(1,N);

for i = 1 : N   
   dp = dx(i,:).*dp + p(i,:);
   end

e   = [1 meshgrid(a+1,ones(N-1,1))'];
K   = [a+1 meshgrid((a+1)^2,ones(N-1,1))'; meshgrid(e,ones(N-1,1))];
dp2 = (dp'*(1./dp)).^2;
w = 1./sum(dp2.*K,2);

% ------------- differentiation arrays ------------------ %

function [A,B] = lDM(x,N);

p  = [ones(N,1) zeros(N,3)];
xt = x(:,ones(2));

for j = 1 : N
   p = (xt-x(j)).*p + [zeros(N,1) p(:,1) 2*p(:,2) 3*p(:,3)];
   end

dx = x(:,ones(1,N)) - x(:,ones(1,N))' + eye(N)*2;
[p1 p2] = meshgrid(p(:,2));
A = (triu(p2,1)+tril(p2,-1)+diag(p(:,3)))./p1./dx;
B = 2*A.*(ndgrid(diag(A))-1./dx);
B = triu(B,1) + tril(B,-1) + diag(p(:,4)./p(:,2)/3);

% ------------- roots of a Jacobi polynomial ------------------ %

function [r] = rJ(b,N);
                      
x = (cos(pi*[N*3-1:-1:0]'/(N*3-1))+1)/2;      
r = rdf( x,dJs(x,b,N,1) );                     
update = r;
j = [1:length(r)]';

while ~isempty(j)
   update(j) = -dJs(r(j),b,N);
   r(j)      = r(j) + update(j);
   j         = find(abs(update)>sqrt(eps));
   end