lobatto.cpp

The Spectral Toolkit is a C++ spectral transform library written by Rodney James and Chuck Panaccion

//
// spectral toolkit 
// copyright (c) 2005 university corporation for atmospheric research
// licensed under the gnu general public license
//

#include "lobatto.h"

using namespace spectral;

/// Constructor for Gauss-Lobatto quadrature.  
/// \param N total length of stencil on \f$ [-1,1] \f$

lobatto::lobatto(int N) : gauss(N)
{
  int m=N-1;

  jpx = new real[m+2];
  djpx= new real[m+2];
    
  real *xj  = new real[m+2];
  real *glw = new real[m+2];
  real *dxx = new real[m+2];
  real *xjr = new real[m+2];
  real *xpt = new real[m+2];
  real *dpt = new real[m+2];
  real *xpt1= new real[m+2];
  real *dpt1= new real[m+2];
    
  real alpha, beta;
  real rr, sm,  djp, xjp;    
  real dth,dd, aa,bb, cd,cs,ss,sd,tmp;     
    
  const int itn = 25;         // Max iteration number
  alpha = 0.0;          
  beta = 0.0;          
    
  xj[0] = 1.0;
  xj[m] = -1.0;  
    
  int nh = (m+1)/2;
    
  jacobi(m+1,alpha,beta, 1.0,xpt, dpt);
  jacobi(m+1,alpha,beta,-1.0,xpt1,dpt1);
    
  dd =  xpt[m]*xpt1[m-1] - xpt1[m]*xpt[m-1];
  aa = (xpt1[m+1]*xpt[m-1] - xpt1[m-1]*xpt[m+1] )/dd;
  bb = (xpt1[m]*xpt[m+1] - xpt1[m+1]*xpt[m] )/dd;
    
  dth = pi/(real)(2*m+1);
  cd = std::cos(2.0*dth);
  sd = std::sin(2.0*dth);
  cs = std::cos(dth);
  ss = std::sin(dth);
    
  for(int k=1;k<nh;k++)
    {
      rr = cs;    
      real dl = 1.0; 
      int iteration=0;
      dl=2.0*eps;
      while((std::abs(dl) > eps)&&(iteration<itn))
        {
	  jacobi(m+1,alpha,beta,rr,xpt,dpt);
	  xjp =  xpt[m+1]+aa* xpt[m]+bb* xpt[m-1];
	  djp =  dpt[m+1]+aa* dpt[m]+bb* dpt[m-1];
	  sm = 0.0;  
	  for(int j=0;j<k;j++)
	    sm = sm + 1.0 /(rr - xj[j]);
	  dl = -(xjp /(djp - xjp*sm));  
	  rr = rr + dl; 
	  iteration++;
        } 
      xj[k] = rr;
      tmp = cs*cd-ss*sd;
      ss = cs*sd+ss*cd;
      cs = tmp;
    }
  for(int k=1;k<=nh;k++)
    xj[m-k] = -xj[k];
  if(m%2==0) 
    xj[nh] = 0.0;   
  for(int k=0;k<m+1;k++)
    {
      rr = xj[k];
      jacobi(m,alpha,beta,rr,xpt,dpt);
      xjp = xpt[m];
      glw[k] = 2.0 / ((real)(m*(m+1)) * xjp*xjp);
      dxx[m-k] = xjp;
    }
  for(int j = 0;j<m+1;j++)
    {
      xjr[m-j] = xj[j];
      points[m-j] = xj[j];
      weights[m-j] = glw[j];
    }
  delete [] xj;
  delete [] glw;
  delete [] dxx;
  delete [] xjr;
  delete [] xpt;
  delete [] dpt;
  delete [] xpt1;
  delete [] dpt1;
}

/// Destructor for Gauss-Lobatto quadrature class.

lobatto::~lobatto()
{
  delete [] jpx;
  delete [] djpx;
}    

/// Internal routine.

void lobatto::jacobi(int nox,real alpha,real beta,real x,real *px,real *dpx)
{
  real k1,dk, ab, ac1,ac2,ac3,acx;
  
  for(int k = 0;k<nox+1;k++)
    {
      px[k] = 0.0;    
      dpx[k] = 0.0;   
    }
  switch(nox)
    {
    case 0:
      {
	px[nox] = 1.0;
	dpx[nox] = 0.0;
	break;
      }
    case 1:
      {
	px[0] = 1.0;
	px[1] =  (1.0 + alpha)*x;
	dpx[0] = 0.0; 
	dpx[1] = 1.0 + alpha;
	break;
      }
    default:
      {
	ab = alpha + beta; 
	jpx[0] = 1.0; 
	jpx[1] = (1.0 + alpha)*x;
	djpx[0] = 0.0; 
	djpx[1] = (1.0 + alpha);
	for(int k = 1;k< nox;k++)
	  {
	    dk = (real)k;
	    k1 = (real)(k+1);
	    ac1 =  2.0*k1* (dk+ ab +1.0) * (2.0*dk + ab);
	    ac2 =  (2.0*k1+ ab) * (2.0*dk +ab + 1.0) * (2.0*dk + ab);
	    acx =  (2.0*dk + ab + 1.0) * (alpha*alpha - beta*beta) + x *ac2; 
	    ac3 =  2.0*(dk + alpha) * (dk + beta) * (2.0*k + ab + 2.0);
	    jpx[k+1] = (acx*jpx[k] - ac3*jpx[k-1]) / ac1;
	    djpx[k+1] = (acx*djpx[k] + ac2*jpx[k] - ac3*djpx[k-1]) / ac1;
	  }
      }
    }
  for(int k = 0;k< nox+1;k++)              
    {
      px[k] = jpx[k];
      dpx[k] = djpx[k];
    }
}