Orthogonal polynomials in gsl/poly

Orthogonal polynomials in gsl/poly

[Richard Mathar]

[-- Attachment #1: Type: text/plain, Size: 661 bytes --]

It would be useful, given the general importance of orthogonal polynomials
in applications, to have some shortcuts to evaluate the most important
ones in "by name" in  the subdirectory poly, instead of defining them by
a list of values that are interpolated.
A proposal/implementation for the Hermite, Laguerre and Chebyshev polynomials
is attached. I can implement more of these if this finds support.

-- 
---
Richard J Mathar                Tel (+31) (0) 71 527 8459
Sterrewacht Universiteit Leiden Fax (+31) (0) 71 527 5819
Postbus 9513
The Netherlands                 URL http://www.strw.leidenuniv.nl/~mathar
office: Niels Bohrweg 2, 2333 CA Leiden, #456

[-- Attachment #2: gsl_poly_orth.h --]
[-- Type: text/plain, Size: 1354 bytes --]

/* poly/gsl_poly_orth.h
*
* 2006-09-09 Richard J. Mathar, Leiden Observatory
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or (at
* your option) any later version.
* 
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
* General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

#ifndef __GSL_POLY_ORTH_H__
#define __GSL_POLY_ORTH_H__


#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

#include <stdlib.h>

int gsl_poly_dd_init_orthH (double dd[], double xa[], const size_t size) ;
int gsl_poly_dd_init_orthL (double dd[], double xa[], const size_t size) ;
int gsl_poly_dd_init_orthT (double dd[], double xa[], const size_t size) ;

__END_DECLS

#endif /* __GSL_POLY_ORTH_H__ */

[-- Attachment #3: gsl_poly_orth.c --]
[-- Type: text/plain, Size: 4628 bytes --]

/* gsl/poly/gsl_poly_orth.c
 *
 * Copyright (C) 2006 Richard J. Mathar
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */


#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_pow_int.h>

/**
* @param[out] dd divided differences for Hermite Polynomial H(n,x).
*   Array must have been allocated with n+1 elements on input and is modified
*   on output.
* @param[out] xa abscissa values for the dd array
*   Array must have been allocated with n+1 elements on input and is modified
*   on output.
* @param[in] size One more than the index n of the polynomial; also the
*        length of the arrays dd and xa.
* A call of the function prepares the arrays dd and xa for calls
* with gsl_poly_dd_eval() or gsl_poly_dd_taylor() .
* @since 2006-09-09
* @author Richard J. Mathar
*/
int gsl_poly_dd_init_orthH (double dd[], double xa[], const size_t size)
{
	int l ;
	const int n=size-1 ;
	double fact= gsl_pow_int(2.,n) ;
	/* We use the representation 25.1.10 of the divided differences,
	* evaluating the derivatives of the Hermite polynomial at x=0.
	*/
	for(l=0 ; l < size ; l++)
		xa[l]=0. ;

	/* Since the polynomial has even parity for even n and odd parity for
	* odd n, the l'th derivative at x=0 is zero for odd n-l.
	*/
	for(l= n-1; l >= 0 ; l -= 2)
		dd[l]=0. ;

	/* the l'th derivative is generally computed by eq 22.3.10 of
	* Abramowitz & Stegun
	*/
	for(l= n; l >= 0 ; l -= 2)
	{
		dd[l]= fact ;
		fact *= -0.5*l*(l-1)/(double)(n-l+2) ;
	}

  	return GSL_SUCCESS;
}

/**
* @param[out] dd divided differences for Laguerre Polynomial L(n,x).
*   Array must have been allocated with n+1 elements on input and is modified
*   on output.
* @param[out] xa abscissa values for the dd array
*   Array must have been allocated with n+1 elements on input and is modified
*   on output.
* @param[in] size One more than the index n of the polynomial; also the
*        length of the arrays dd and xa.
* After a call of the function, the arrays dd and xa can be used for calls
* with gsl_poly_dd_eval() or gsl_poly_dd_taylor() .
* @since 2006-09-09
* @author Richard J. Mathar
*/
int gsl_poly_dd_init_orthL (double dd[], double xa[], const size_t size)
{
	int l ;
	const int n=size-1 ;
	double fact= 1.0 ;
	/* We use the representation 25.1.10 of the divided differences,
	* evaluating the derivatives of the Laguerre polynomial at x=0.
	*/
	for(l=0 ; l < size ; l++)
		xa[l]=0. ;

	/* the l'th derivative is generally computed by eq 22.3.9 of
	* Abramowitz & Stegun for the case alpha=0.
	*/
	for(l= 0; l <= n ; l ++)
	{
		dd[l]= fact ;
		fact *= (l-n)/gsl_pow_2(1+l) ;
	}

  	return GSL_SUCCESS;
}

/**
* @param[out] dd divided differences for Chebyshev Polynomial of the 1st kind T(n,x).
*   Array must have been allocated with n+1 elements on input and is modified
*   on output.
* @param[out] xa abscissa values for the dd array
*   Array must have been allocated with n+1 elements on input and is modified
*   on output.
* @param[in] size One more than the index n of the polynomial; also the
*        length of the arrays dd and xa.
* After a call of the function, the arrays dd and xa can be used for calls
* with gsl_poly_dd_eval() or gsl_poly_dd_taylor() .
* @since 2006-09-09
* @author Richard J. Mathar
*/
int gsl_poly_dd_init_orthT (double dd[], double xa[], const size_t size)
{
	int l ;
	const int n=size-1 ;
	double fact= gsl_pow_int(2.,n-1) ;
	/* We use the representation 25.1.10 of the divided differences,
	* evaluating the derivatives of the Chebyshev polynomial at x=0.
	*/
	for(l=0 ; l < size ; l++)
		xa[l]=0. ;

	/* Since the polynomial has even parity for even n and odd parity for
	* odd n, the l'th derivative at x=0 is zero for odd n-l.
	*/
	for(l= n-1; l >= 0 ; l -= 2)
		dd[l]=0. ;

	/* the l'th derivative is generally computed by eq 22.3.6 of
	* Abramowitz & Stegun, here only at x=0.
	*/
	for(l= n; l >= 0 ; l -= 2)
	{
		dd[l]= fact ;
		fact *= l*(1-l)/(double)((n-l+2)*(n+l-2)) ;
	}

  	return GSL_SUCCESS;
}

Re: Orthogonal polynomials in gsl/poly

[Brian Gough]

At Sat, 9 Sep 2006 18:02:58 +0200,
Richard Mathar wrote:
> It would be useful, given the general importance of orthogonal polynomials
> in applications, to have some shortcuts to evaluate the most important
> ones in "by name" in  the subdirectory poly, instead of defining them by
> a list of values that are interpolated.
> A proposal/implementation for the Hermite, Laguerre and Chebyshev polynomials
> is attached. I can implement more of these if this finds support.

Thanks for the email.  I agree it would be useful to access the
different polynomials by name, there are a few in specfunc (based on
explicit representations for small n and recurrences for general n I
think) but more would be useful.  A first question would be regarding
the relative merits of different ways of computing them, e.g stability
& speed, as I haven't looked at this problem so I have no idea what
the best way is.

-- 
Brian Gough