Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

Function.cc

//
// some particular functions on [0,1]^D
//

# include <stdlib.h>
# include <math.h>

# include "Function.hpp"

double Function::GaussQuadrature(double *XA, double *XE, int D, int order) {
 double x[100] ;
 if ((order&1)!=0)  { std::cout<< "Function::GaussQuadrature: order must be even\n"; exit(-1) ;}
 return GQ(x,0,XA,XE,D , order) ;
} 


double Function::GQ(double *x, int i, double *XA, double *XE, int D, int p) {
 if (i==D) return Eval(x) ;
 
 int j;
 double w[100],n[100] ;
 switch (p) {
   case 4: {
     w[0]= 3.47854845137454e-01, n[0]=-8.61136311594052e-01 ;
     w[1]= 6.52145154862546e-01, n[1]=-3.39981043584856e-01 ;
    }
   break ;
   case 6: {
     w[0]= 1.71324492379171e-01, n[0]= -9.32469514203152e-01;
     w[1]= 3.60761573048139e-01, n[1]= -6.61209386466264e-01;
     w[2]= 4.67913934572691e-01, n[2]= -2.38619186083197e-01;
    }
   break ; 
   case 8: {
     w[0]= 1.01228536290376e-01, n[0]= -9.60289856497536e-01;
     w[1]= 2.22381034453375e-01, n[1]= -7.96666477413627e-01;
     w[2]= 3.13706645877887e-01, n[2]= -5.25532409916329e-01;
     w[3]= 3.62683783378362e-01, n[3]= -1.83434642495649e-01; 
   }
   break ;
   case 10: {
     w[0]= 6.66713443086885e-02, n[0]= -9.73906528517172e-01;
     w[1]= 1.49451349150580e-01, n[1]= -8.65063366688985e-01;
     w[2]= 2.19086362515982e-01, n[2]= -6.79409568299024e-01;
     w[3]= 2.69266719309997e-01, n[3]= -4.33395394129247e-01;
     w[4]= 2.95524224714753e-01, n[4]= -1.48874338981631e-01;
   }
   break ;
   case 12: {
     w[0]= 4.71753363865110e-02, n[0]= -9.81560634246720e-01;
     w[1]= 1.06939325995318e-01, n[1]= -9.04117256370475e-01;
     w[2]= 1.60078328543346e-01, n[2]= -7.69902674194305e-01;
     w[3]= 2.03167426723066e-01, n[3]= -5.87317954286618e-01;
     w[4]= 2.33492536538355e-01, n[4]= -3.67831498998180e-01;
     w[5]= 2.49147045813403e-01, n[5]= -1.25233408511469e-01;
   }
   break ;
   case 14: {
    w[0]=  3.51194603317506e-02; n[0]= -9.86283808696813e-01;
    w[1]=  8.01580871597603e-02; n[1]= -9.28434883663574e-01;
    w[2]=  1.21518570687903e-01; n[2]= -8.27201315069765e-01;
    w[3]=  1.57203167158194e-01; n[3]= -6.87292904811685e-01;
    w[4]=  1.85538397477938e-01; n[4]= -5.15248636358154e-01;
    w[5]=  2.05198463721296e-01; n[5]= -3.19112368927890e-01;
    w[6]=  2.15263853463158e-01; n[6]= -1.08054948707344e-01;
   }
   break ;
   case 16: {
    w[0]=  2.71524594117552e-02; n[0]= -9.89400934991649e-01;
    w[1]=  6.22535239386477e-02; n[1]= -9.44575023073233e-01;
    w[2]=  9.51585116824926e-02; n[2]= -8.65631202387832e-01;
    w[3]=  1.24628971255534e-01; n[3]= -7.55404408355004e-01;
    w[4]=  1.49595988816577e-01; n[4]= -6.17876244402644e-01;
    w[5]=  1.69156519395003e-01; n[5]= -4.58016777657228e-01;
    w[6]=  1.82603415044924e-01; n[6]= -2.81603550779259e-01;
    w[7]=  1.89450610455069e-01; n[7]= -9.50125098376373e-02;
   }
   break ;
 }

 for (j=p/2; j<p; j++) {w[j]=w[p-1-j]; n[j]=-n[p-1-j] ;}

 double f=0.0,DX2=(XE[i]-XA[i])/2 ;
 for (j=0; j<p; j++) {
   x[i]=XA[i]+DX2*(n[j]+1) ;
   f +=w[j]*DX2*GQ(x,i+1,XA,XE,D,p) ;
 }
 return f ;
}

---