00001 #include <math.h>
00002 #include "tbl.h"
00003 #include <assert.h>
00004
00005 double* jacobi(int , double) ;
00006
00007 double* pointsgausslobatto(int n) {
00008
00009 int nmax = 10 ;
00010 int ndiv = (n+1)*(n+1)+5 ;
00011 double pas = double(2)/double(ndiv-1) ;
00012 double eps = 2.5E-12 ;
00013
00014 Tbl subdiv(ndiv) ;
00015 subdiv.set_etat_qcq();
00016 Tbl cs(n-1,2) ;
00017 cs.set_etat_qcq() ;
00018 double* xx = new double[nmax] ;
00019 double* yy ;
00020 double* zz ;
00021 int i,k,l ;
00022
00023 for (i = 0 ; i < ndiv ; i++) {
00024 subdiv.set(i) = double(-1) + pas * double(i) ;
00025 }
00026
00027 double a = (2*double(n)+3)/double((n+1)*(n+1)) ;
00028 double b = - 1 - a ;
00029
00030 int j=0;
00031
00032 for (i = 1 ; i < ndiv-3 ; i++) {
00033 yy = jacobi(n+1 , subdiv(i)) ;
00034 zz = jacobi(n+1 , subdiv(i+1)) ;
00035 double omega1 = yy[n+1] + a * yy[n] + b * yy[n-1] ;
00036 double omega2 = zz[n+1] + a * zz[n] + b * zz[n-1] ;
00037 if (omega1*omega2 <= 0) {
00038 cs.set(j,0) = subdiv(i) ;
00039 cs.set(j,1) = subdiv(i+1) ;
00040 j++;
00041 }
00042 delete [] yy ;
00043 delete [] zz ;
00044 }
00045
00046
00047
00048 double* pointsgl = new double[j+2];
00049 assert(j==n-1) ;
00050 pointsgl[0] = -1;
00051
00052 for (l = 0 ; l < j ; l++) {
00053 xx[0] = cs(l,0) ;
00054 xx[1] = cs(l,1) ;
00055 for (k = 2 ; k < nmax ; k++) {
00056 yy = jacobi(n+1 , xx[k-2]) ;
00057 zz = jacobi(n+1 , xx[k-1]) ;
00058 double omega1 = yy[n+1] + a * yy[n] + b * yy[n-1] ;
00059 double omega2 = zz[n+1] + a * zz[n] + b * zz[n-1] ;
00060 if (( fabs(xx[k-2]-xx[k-1])>=eps )&&(fabs(omega2)>=eps )) {
00061 xx[k]=(xx[k-2]*omega2 -xx[k-1]*omega1)/(omega2-omega1) ;
00062 }
00063 else {
00064 xx[k]=xx[k-1];
00065 }
00066 delete [] yy ;
00067 delete [] zz ;
00068 }
00069 pointsgl[l+1] = xx[nmax-1] ;
00070 }
00071 delete [] xx ;
00072
00073 pointsgl[n] = 1 ;
00074
00075 return pointsgl ;
00076 }
