fejer2.cpp

/*
 * Copyright (c) 2005, Andrew Fernandes (andrew@fernandes.org);
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 
 * - Redistributions of source code must retain the above copyright
 * notice, this list of conditions and the following disclaimer.
 * 
 * - Redistributions in binary form must reproduce the above copyright
 * notice, this list of conditions and the following disclaimer in the
 * documentation and/or other materials provided with the distribution.
 * 
 * - Neither the name of the North Carolina State University nor the
 * names of its contributors may be used to endorse or promote products
 * derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */
#include "gaussqr.h"
#include "linalg.hpp"
#include "fftpack.h"

#include <new>
using std::nothrow;

#include <cmath>
using std::cos;

//****************************************************************************80

gaussqr_result fejer2_abscissae ( const integer_t n , real_t *z , real_t *q )

//****************************************************************************80
//
//  Purpose:
//
//    FEJER2_ABSCISSAE computes points and weights of a Fejer Type 2 rule.
//
//  Discussion:
//
//    The rule is defined for the interval [-1,+1].
//
//    Weights q and nodes z for the n-point Fejer type-2 quadrature rules,
//    calculated via the inverse (fast) real_t discrete fourier transform.
//
//    The knots are the same knots as the Clenshaw-Curtis knots, implying that
//    by doubling the number of evalulation points, you can re-use the function
//    evaluations at the old points, since the new set of knots will contain the
//    old set of knots. In other words, the knots for an n-point fejer2 rule
//    will be contained in the (2*n+1)-point fejer2 rule.
// 
//    Although capable of computing knots and weights for arbitrary 'n', the
//    computation will be much faster and efficient if (n+1) has factors only
//    from the set {2,3,4,5}.
//    
//    Therefore one recommended sequence of n is {3,7,15,31,63,127,255,511,1023,...}.
//    Three (or five) times this set would also be acceptable.
//
//  Modified:
//
//    29 April 2013
//
//  Author:
//
//    Andrew Fernandes
//
//  Parameters:
//
//    Input, int N, the order of the rule.
//    1 <= N.
//
//    Output, real Z[N], Q[N], the points and weights of the rule.  
//
{
  if ( n == 1 )
  {
    z[0] = 0.0;
    q[0] = 2.0;
    return gaussqr_success;
  }

  if ( n == 2 )
  {
    z[0] = - 0.5;
    z[1] = + 0.5;
    q[0] = 1.0;
    q[1] = 1.0;
    return gaussqr_success;
  }

  if ( z == 0 )
  {
    return gaussqr_illegal_argument;
  }
  if ( q == 0 )
  {
    return gaussqr_illegal_argument;
  }
    
  gaussqr_result rv = gaussqr_success;

  integer_t n_fft = n + 1;
  const real_t pi = 3.1415926535897932384626433832795028841971693993751058209749445923;
//
//  Allocate memory
//
  real_t *v = new(nothrow) real_t[n_fft];
  real_t *work = new(nothrow) real_t[2*n_fft];
  integer_t *ifac = new(nothrow) integer_t[sizeof(integer_t)*8];

  if ( v == 0 || work == 0 || ifac == 0 ) 
  {
    rv = gaussqr_memory_allocation_error;
    goto done;
  }
//
// calculate the knots
//
  for ( integer_t k = 1; k < n_fft; k++ ) 
  {
    z[k-1] = -cos((k*pi)/n_fft);
  }
//
// calculate the transformed weights
//
  v[0] = 2.0;
  for ( integer_t k = 1; k < n_fft; k++ )
  {
    v[k] = 0.0;
  }
  for ( integer_t k = 1; k < n_fft/2; k++ ) 
  {
    v[2*k-1] = 2.0/(1.0-4.0*(k*k));
  }
  v[n_fft-1-(n_fft&1)] = (n_fft-3.0) / (2.0*(n_fft/2)-1.0) - 1.0;
//
//  Initialize the FFT package.
//
  rffti ( &n_fft, work, ifac );
//
// take the fourier transform, and scale
//
  rfftb ( &n_fft, v, work, ifac );

  for ( integer_t i = 0; i < n_fft; i++ )
  {
    v[i] /= n_fft;
  }
//
// transfer u to q, removing first (zero) element
//
  for ( integer_t i = 1; i < n_fft; i++ ) 
  {
    q[i-1] = v[i];
  }

done:
    delete[](ifac);
    delete[](work);
    delete[](v);

    return(rv);
}


gaussqr_result map_fejer2_domain( const real_t a, const real_t b , const domain_type type , 
  const integer_t n , const real_t *x , real_t *y , real_t *dy )
//
//    Given an interval [a,b] and a domain_type, this subroutine maps the interval [-1,1] to [a,b],
//    an interval that may be closed, left-infinite, right-infinite, or totally infinite.
//
//    The input points are x, and should be enclosed in the open interval (-1:1). 
//    This restriction is NOT checked for.
// 
//    On output, the points y are the mapped points to the new interval,
//    and the dy array gives the derivative of the transformation.
//
//    Note that the change of integration variables is 
//    $\int_{-1}^{+1} f(x) \, dx = \int_{a}^{b} f(y(x)) dy \, dx$, and this is
//    used in the quadrature approximations to integration.
//
{
  if ( a >= b || n < 3 || x == 0 || y == 0 || dy == 0 )
        return(gaussqr_illegal_argument);
    
    switch(type) {
        
        case domain_finite :
            for ( integer_t i = 0; i < n; i++ ) {
                const real_t bpa = b + a;
                const real_t bma = b - a;
                y[i] = 0.5 * ( x[i]*bma + bpa );
                dy[i] = 0.5 * bma;
            }
            break;
    
        case domain_left_infinite :
            for ( integer_t i = 0; i < n; i++ ) {
                const real_t x_m1 = (1.0-x[i]);
                const real_t x_p1 = (1.0+x[i]);
                y[i] = x_m1/x_p1 + b;
                dy[i] = 2.0/(x_p1*x_p1);
            }
            break;
            
        case domain_right_infinite :
            for ( integer_t i = 0; i < n; i++ ) {
                const real_t x_m1 = (1.0-x[i]);
                const real_t x_p1 = (1.0+x[i]);
                y[i] = x_p1/x_m1 + a;
                dy[i] = 2.0/(x_m1*x_m1);
            }
            break;
            
        case domain_infinite :
            for ( integer_t i = 0; i < n; i++ ) {
                const real_t x_m1 = (1.0-x[i]);
                const real_t x_p1 = (1.0+x[i]);
                y[i] = (2.0*x[i])/(x_m1*x_p1);
                dy[i] = (2.0*(1.0+(x[i]*x[i])))/(x_m1*x_m1*x_p1*x_p1);
            }
            break;
            
        default :
            return(gaussqr_illegal_argument);
        
    }
    
    return(gaussqr_success);
}