sandbox/easystab/dif1D.m

    Function dif1D

    function [d,dd,w,s]=dif1D(type,s0,L,N,opt);
    

    This function creates a mesh s, differentiation matrices of first and second order d and dd, as well as the weight w.

    On output: s is a column-vector containing the mesh points, starting from the point at the left. d and dd are the derivation matrices (quare matrices, in full or sparse storing format) w is a line-vector containing the weights.

    Last parameter ‘opt’ is optional. For Finite-Difference it is interpreted as the number of points of the stencil (default :3) and for Stretched Chebyschev it is interpreted as the parameter b controlling the stretching.

    switch type
    

    Finite differences (optimized case)

    For a domain [s0,s0+L]

    Last parameter ‘pts’ is optional, if absent it is set to 3 (relevant value in almost all cases)

       case 'fd' 
           
      if(nargin==4) 
        pts = 3;
      else
       pts = opt;
      end
      scale=L/2;
      [s,d] = fddif(N,1,pts); 
      [s,dd] = fddif(N,2,pts); 
      s=s*scale; 
      d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
      w=([diff(s'),0]+[0,diff(s')])/2;    
    

    Finite differences (basic implementation, equivalent to the previous)

     case 'fds' 
        
    
      dx = L/(N-1); % grid spacing
      s = [  s0 : L/(N-1) : s0+L ]';  % creates an equispaced mesh with N point on the interval
      for i=2:N-1 % centered three-point formulas 
        d(i,i) = 0;
        d(i,i+1) =  1/(2*dx); 
        d(i,i-1) = -1/(2*dx);
        dd(i,i) = -2/dx^2;
        dd(i,i+1) =  1/(dx^2); 
        dd(i,i-1) =  1/(dx^2);        
      end
     % uncentred formulas for first and last grid points 
       d(1,1) = -3/(2*dx); d(1,2) = 4/(2*dx) ; d(1,3) = -1/(2*dx);
       d(N,N-2) =  1/(2*dx); d(N,N-1) = -4/(2*dx); d(N,N)   =  3/(2*dx);
       dd(1,1) = 2/dx^2; dd(1,2) = -5/dx^2; dd(1,3) = 4/dx^2; dd(1,4) = -1/dx^2; 
       dd(N,N) = 2/dx^2; dd(N,N-1) = -5/dx^2; dd(N,N-2) = 4/dx^2; dd(N,N-3) = -1/dx^2; 
      % "weight" to compute integrals using the trapeze rule
      for i=2:N-1
        w(i) = dx;
      end
      w(1) = dx/2; w(N) = dx/2;
    

    Finite differences on a periodic domain

    Domain is defined as [s0,s0+L]

     case 'fp' % periodic finite differences 
      scale=L/2;
      [s,d] = fdper(N,1,opt); 
      [s,dd] = fdper(N,2,opt);
      s=s*scale; 
      d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
      w=ones(1,N)*(s(2)-s(1)); 
    

    Chebyshev collocation

    on domain [s0,s0+L]

     case 'cheb' % chebychev 
      scale=-L/2;
      [s,DM] = chebdif(N,2);
      d=DM(:,:,1);  
      dd=DM(:,:,2);
      s=s*scale; 
      d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
      w=L*clencurt(N)/2; 
    

    Fourier collocation

    (requires a periodic domain [s0,s0+L])

      case 'fou' % Fourier 
      scale=L/(2*pi); 
      [s, d] = fourdif(N, 1);
      [s, dd] = fourdif(N, 2);
      s=s*scale; 
      w=ones(1,N)*(s(2)-s(1));
      d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
    

    Hermite collocation on infinite domain.

    Domain is centred on s0, L controls the density

    (increase L to cluster the points towards s0)

      case 'her' % hermite 
      scale=1;
      [s,DM] = herdif(N,2,L);
      d=DM(:,:,1);  
      dd=DM(:,:,2);
      s=s*scale; 
      d=d/scale; dd=dd/scale^2; s=s+s0;
      w=zeros(1,N); % weight is not implemented yet !!!
    

    Stretched Chebyshev ; infinite ; exponential

    • s0 is the center of the domain.
    • L controls the clustering of points around the origin : about half points will be localised in [s0-L,s0+L]
    • optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
      case 'chebInfExp' % chebychev with coordinate stretching to treat large domains
          
     
       if(nargin==4) 
          b = 0.999;
       else
          b = opt;
       end
       
       [dxi,ddxi,wxi,xi]=dif1D('cheb',-b,2*b,N); % Standart Chebyshev for xi \in [-b,b]
       
      s = L*atanh(xi)+s0; % coordinate stretching (exponential)
      G(:, 1) =                          - (xi.^2 - 1) /L;      % G1 = ds / dxi
      G(:, 2) =                   2 *xi .* (xi.^2 - 1) /L^2;    % G2 = d^2 s /d xi^2
      G(:, 3) =    - 2 * (3 * xi.^2 - 1) .* (xi.^2 - 1) / L^3;
      G(:, 4) = 8 * xi .* (3 * xi.^2 - 2) .* (xi.^2 - 1) / L^4;
      
      d = G(:,1).*dxi;
      dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
      w = wxi./G(:,1)';
    

    Stretched Chebyshev ; infinite ; algebraic

    • s0 is the center of the domain.
    • L controls the clustering of points around the origin : about half points will be localised in [s0-L,s0+L]
    • optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
      case 'chebInfAlg' % chebychev with coordinate stretching to treat large domains
          
      if(nargin==4) 
          b = 0.999;
       else
          b = opt;
      end
      
      [dxi,ddxi,wxi,xi]=dif1D('cheb',-b,2*b,N); % Standart Chebyshev for xi \in [-b,b]
      
      s = L * xi ./ sqrt(1 - xi.^2); 
      K = sqrt(s.^2 + L^2);
      G(:, 1) =                                    L^2 ./ K.^3;
      G(:, 2) =                          - 3 * s * L^2 ./ K.^5;
      G(:, 3) =             3 * (4 * s.^2 - L^2) * L^2 ./ K.^7;
      G(:, 4) = - 15 * s .* (4 * s.^2 - 3 * L^2) * L^2 ./ K.^9;
     
      d = G(:,1).*dxi;
      dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
      w = wxi./G(:,1)';
    

    Stretched Chebyshev ; semi-infinite ; exponential

    • s0 is the left bound of the domain.
    • L controls the clustering of points around the origin : about half points will be localised in [s0,s0+L]
    • optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
     case 'chebSemiInfExp' % chebychev with coordinate stretching to treat large domains
     
     if(nargin==4) 
          b = 0.999;
       else
          b = opt;
     end
         
      [dxi,ddxi,wxi,xi]=dif1D('cheb',0,b,N); % Standart Chebyshev for xi \in [0,b]
      
      s = L*atanh(xi)+s0; % coordinate stretching (exponential)
      G(:, 1) =                          - (xi.^2 - 1) /L;      % G1 = ds / dxi
      G(:, 2) =                   2 * xi .* (xi.^2 - 1) /L^2;    % G2 = d^2 s /d xi^2
      d = G(:,1).*dxi;
      dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
      w = wxi./G(:,1)';
    

    Stretched Chebyshev ; semi-infinite ; algebraic

    • s0 is the left bound of the domain.
    • L controls the clustering of points around the origin : about half points will be localised in [s0,s0+L]
    • optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
      case 'chebSemiInfAlg' % chebychev with coordinate stretching to treat large domains
          
      if(nargin==4)  
          b = 0.999;  
      else
          b = opt; 
      end
       
      [dxi,ddxi,wxi,xi]=dif1D('cheb',0,b,N); % Standart Chebyshev for xi \in [0,1]
    
      s = L * xi ./ sqrt(1 - xi.^2); 
      K = sqrt(s.^2 + L^2);
      G(:, 1) =                                    L^2 ./ K.^3;
      G(:, 2) =                          - 3 * s .* L^2 ./ K.^5;
     
      d = G(:,1).*dxi;
      dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
      w = wxi./G(:,1)';
      
    
    end