00001 /* 00002 * Definition of class Chebyshev_poly 00003 */ 00004 00005 /* 00006 * Copyright (c) 2005 Eric Gourgoulhon 00007 * 00008 * This file is part of LORENE. 00009 * 00010 * LORENE is free software; you can redistribute it and/or modify 00011 * it under the terms of the GNU General Public License as published by 00012 * the Free Software Foundation; either version 2 of the License, or 00013 * (at your option) any later version. 00014 * 00015 * LORENE is distributed in the hope that it will be useful, 00016 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00017 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00018 * GNU General Public License for more details. 00019 * 00020 * You should have received a copy of the GNU General Public License 00021 * along with LORENE; if not, write to the Free Software 00022 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00023 * 00024 */ 00025 00026 /* 00027 * $Id: chebyshev_poly.C,v 1.2 2005/11/14 14:12:09 e_gourgoulhon Exp $ 00028 * $Log: chebyshev_poly.C,v $ 00029 * Revision 1.2 2005/11/14 14:12:09 e_gourgoulhon 00030 * Added include <assert.h> 00031 * 00032 * Revision 1.1 2005/11/14 01:56:58 e_gourgoulhon 00033 * First version 00034 * 00035 * 00036 * $Header: /cvsroot/Lorene/School05/Monday/chebyshev_poly.C,v 1.2 2005/11/14 14:12:09 e_gourgoulhon Exp $ 00037 * 00038 */ 00039 00040 00041 #include <iostream> 00042 00043 using namespace std ; 00044 00045 #include <stdlib.h> 00046 #include <math.h> 00047 #include <assert.h> 00048 00049 #include "ortho_poly.h" 00050 #include "grid.h" 00051 00052 void coef_cheb_fft(int np, const double* ff, double* cf) ; 00053 00054 //----------------------// 00055 // Standard constructor // 00056 //----------------------// 00057 00058 Chebyshev_poly::Chebyshev_poly(int ni) : Ortho_poly(ni) {} 00059 00060 //--------------------// 00061 // Copy constructor // 00062 //--------------------// 00063 00064 Chebyshev_poly::Chebyshev_poly(const Chebyshev_poly& pp) : Ortho_poly(pp) {} 00065 00066 00067 //--------------// 00068 // Destructor // 00069 //--------------// 00070 00071 Chebyshev_poly::~Chebyshev_poly() {} 00072 00073 00074 //--------------// 00075 // Assignment // 00076 //--------------// 00077 00078 void Chebyshev_poly::operator=(const Chebyshev_poly& ) { 00079 00080 cerr << "Chebyshev_poly::operator= not implemented !" << endl ; 00081 abort() ; 00082 00083 } 00084 00085 //-----------// 00086 // Display // 00087 //-----------// 00088 00089 ostream& operator<<(ostream& ost, const Chebyshev_poly& pp) { 00090 00091 ost << "Basis of Chebyshev polynomials up to degree " << pp.n() << endl ; 00092 00093 return ost ; 00094 } 00095 00096 //-------------------// 00097 // Weight function // 00098 //-------------------// 00099 00100 double Chebyshev_poly::weight(double x) const { 00101 00102 return 1. / sqrt( 1. - x*x) ; 00103 00104 } 00105 00106 00107 //-------------------------------------------// 00108 // Evaluation of the Chebyshev polynomials // 00109 //-------------------------------------------// 00110 00111 double Chebyshev_poly::operator()(int i, double x) const { 00112 00113 assert( i >= 0 ) ; 00114 assert( i <= nn ) ; 00115 00116 if (i==0) return 1. ; 00117 00118 if (i==1) return x ; 00119 00120 double tjm2 = 1. ; // value at step j - 2 00121 double tjm1 = x ; // value at step j - 1 00122 00123 for (int j=2; j<=i; j++) { 00124 double tj = 2*x* tjm1 - tjm2 ; 00125 tjm2 = tjm1 ; 00126 tjm1 = tj ; 00127 } 00128 00129 return tjm1 ; 00130 00131 } 00132 00133 //------------------------------------// 00134 // Computation of nodes and weights // 00135 //------------------------------------// 00136 00137 00138 const Grid& Chebyshev_poly::gauss_nodes() const { 00139 00140 if (p_gauss_nodes == 0x0) { // the nodes must initialized 00141 00142 p_gauss_nodes = new Grid_Chebyshev_Gauss(nn+1) ; 00143 } 00144 00145 return *p_gauss_nodes ; 00146 } 00147 00148 00149 double Chebyshev_poly::gauss_weight(int i) const { 00150 00151 if (p_gauss_weights == 0x0) { // the weights must be computed 00152 00153 p_gauss_weights = new double[nn+1] ; 00154 00155 for (int j=0; j<=nn; j++) p_gauss_weights[j] = M_PI / double(nn+1) ; 00156 } 00157 00158 return p_gauss_weights[i] ; 00159 00160 } 00161 00162 00163 double Chebyshev_poly::gauss_gamma(int i) const { 00164 00165 if (p_gauss_gamma == 0x0) { // the weights must be computed 00166 00167 p_gauss_gamma = new double[nn+1] ; 00168 00169 p_gauss_gamma[0] = M_PI ; 00170 00171 for (int j=1; j<=nn; j++) p_gauss_gamma[j] = 0.5 * M_PI ; 00172 } 00173 00174 return p_gauss_gamma[i] ; 00175 00176 } 00177 00178 00179 const Grid& Chebyshev_poly::gauss_lobatto_nodes() const { 00180 00181 if (p_gauss_lobatto_nodes == 0x0) { // the nodes must initialized 00182 00183 p_gauss_lobatto_nodes = new Grid_Chebyshev_GL(nn+1) ; 00184 } 00185 00186 return *p_gauss_lobatto_nodes ; 00187 } 00188 00189 00190 double Chebyshev_poly::gauss_lobatto_weight(int i) const { 00191 00192 if (p_gauss_lobatto_weights == 0x0) { // the weights must be computed 00193 00194 p_gauss_lobatto_weights = new double[nn+1] ; 00195 00196 p_gauss_lobatto_weights[0] = M_PI / double(2*nn) ; 00197 00198 for (int j=1; j<nn; j++) 00199 p_gauss_lobatto_weights[j] = M_PI / double(nn) ; 00200 00201 p_gauss_lobatto_weights[nn] = M_PI / double(2*nn) ; 00202 } 00203 00204 return p_gauss_lobatto_weights[i] ; 00205 00206 } 00207 00208 00209 double Chebyshev_poly::gauss_lobatto_gamma(int i) const { 00210 00211 if (p_gauss_lobatto_gamma == 0x0) { // the weights must be computed 00212 00213 p_gauss_lobatto_gamma = new double[nn+1] ; 00214 00215 p_gauss_lobatto_gamma[0] = M_PI ; 00216 00217 for (int j=1; j<nn; j++) p_gauss_lobatto_gamma[j] = 0.5 * M_PI ; 00218 00219 p_gauss_lobatto_gamma[nn] = M_PI ; 00220 00221 } 00222 00223 return p_gauss_lobatto_gamma[i] ; 00224 00225 } 00226 00227 00228 //---------------------------------------------------// 00229 // Gauss-Lobatto interpolant polynomial via a FFT // 00230 //---------------------------------------------------// 00231 00232 void Chebyshev_poly::coef_interpolant_GL_FFT(double (*f)(double), 00233 double* cf) const { 00234 00235 // Values of the function at the Gauss-Lobatto nodes 00236 00237 double* ff = new double[nn+1] ; 00238 00239 for (int i=0; i<=nn; i++) ff[i] = f( gauss_lobatto_nodes()(i) ) ; 00240 00241 // Chebyshev transform via a FFT 00242 00243 coef_cheb_fft(nn+1, ff, cf) ; 00244 00245 delete [] ff ; 00246 00247 } 00248 00249 00250 //---------------------------------------------// 00251 // Coefficient of the orthogonal projection // 00252 //---------------------------------------------// 00253 00254 void Chebyshev_poly::coef_projection(double (*f)(double), double* cf) const { 00255 00256 // The computation is an approximate one: 00257 // it returns the coefficients of an interpolating polynomial with$ 00258 // a large number of Gauss-Lobatto nodes 00259 00260 int n_large = 128 ; 00261 00262 if (nn > n_large) { 00263 cerr << "Chebyshev_poly::coef_projection : nn > n_large !" 00264 << endl << " nn = " << nn << " n_large = " << n_large << endl ; 00265 abort() ; 00266 } 00267 00268 Chebyshev_poly cheb_large(n_large) ; 00269 00270 double* cf_large = new double[n_large + 1] ; 00271 00272 cheb_large.coef_interpolant_GL(f, cf_large) ; 00273 00274 for (int i=0; i<=nn; i++) cf[i] = cf_large[i] ; 00275 00276 delete [] cf_large ; 00277 00278 } 00279 00280 00281 00282 00283 00284 00285 00286 00287 00288 00289 00290