SparseIntegralError.cc

00001 // This program computes the integral expression f(x) := \int_y K(x,y)u(y) dy using
00002 // optimized approximation spaces for K and u according to
00003 // S.Knapek, F. Koster  
00004 // Integral operators on sparse grids. SIAM J. Num. Anal., 39(5):1794-1809, 2002.
00005 
00006 #include "LevelAdaptive.hpp"
00007 #include "IntegralTestFunctions.hpp"
00008 
00009 int debugRefine ;
00010 int main() {
00011   Wavelets WC("Daubechies",3) ;
00012 
00013   // test functions
00014   IntegralTestFunction<2> FU;
00015   IntegralTestFunction<4> FK;
00016   
00017   // functions which define the level adaptive approximation spaces 
00018   SparseGridFunction<2> SGU ;
00019   SparseGridFunction<4> SGK ;
00020   
00021   // boundary condition identifiers for U and K 
00022   int BCU[2][2],BCK[4][2], i ;
00023   for (i=0;i<2; i++) { BCU[i][0]=-1 , BCU[i][1]=-1 ;}
00024   for (i=0;i<4; i++) { BCK[i][0]=-1 , BCK[i][1]=-1 ;}
00025   
00026   
00027   // level 10 uniform grid for error evaulation
00028   int LL0=10 ;
00029   int a[2]={0}, e[2]={1<<LL0 , 1<<LL0} ;
00030   UniformData<2> M(a,e,&WC),M1(a,e,&WC) ;
00031   M.SetBoundaryConditions(BCU) ;
00032   M1.SetBoundaryConditions(BCU) ;
00033   
00034   double l2er,l2era=1.0,h1er,h1era=1.0 ;
00035   
00036   printf("Output:\n") ;
00037   printf("Error <L> <L_2 error> <ratio L_2 error> | <H_1 error> <ratio H_1 error>\n") ;
00038   
00039   for (int LL=4; LL<=9; LL++) {
00040     
00041     int L[4]={LL,LL,LL,LL}, Level0[4]={WC.Level0 , WC.Level0 , WC.Level0 , WC.Level0} ;
00042     
00043     // sparse grid data structures
00044     LevelAdaptiveGrid<2> AU(Level0,L,&WC)  ;
00045     LevelAdaptiveGrid<4> AK(Level0,L,&WC)  ;
00046     BlendingWeights<2>   BU(&AU) ;
00047     BlendingWeights<4>   BK(&AK) ;
00048     
00049     LevelAdaptiveData<2> U(&AU),F(&AU) ;
00050     LevelAdaptiveData<4> K(&AK) ;
00051     
00052     BlendingData<2>      BSU(&BU) ;
00053     BlendingData<4>      BSK(&BK) ;
00054     
00055     // set/allocate level adaptive grids
00056     AU.Set(&SGU,L[0]+(2-1)*WC.Level0) ;
00057     AK.Set(&SGK,L[0]+(4-1)*WC.Level0) ;
00058     
00059     // generate wavelet representations of K and u
00060     // by means of the blending scheme 
00061     BSU.SetFunction(&FU) ;
00062     BSK.SetFunction(&FK) ;
00063     BSU.SetBoundaryConditions(BCU) ;
00064     BSK.SetBoundaryConditions(BCK) ;
00065     
00066     for (i=0; i<2; i++) { 
00067       BSU.ApplyOp(BCU[i] , &BSU , i , POLYQUAD) ;
00068       BSU.ApplyOp(BCU[i] , &BSU , i , TRANSFORM) ;
00069     }
00070     
00071     for (i=0; i<4; i++) {
00072       BSK.ApplyOp(BCK[i] , &BSK , i , POLYQUAD) ;
00073       BSK.ApplyOp(BCK[i] , &BSK , i , TRANSFORM) ;
00074     }
00075     
00076     U.Blend(&BSU) ;
00077     K.Blend(&BSK) ;
00078     
00079     // evaluation of integral operator
00080     K.Mul(0.5*210*210) ;
00081     F.Multiply(&K,&U)  ;
00082     
00083     // nodal values of numerical result on a uniform grid
00084     F.ToUniform(&M) ;
00085     for (i=0; i<2; i++) {
00086       M.ApplyOp(M.Ext.BC[i] , &M , i , INVERSETRANSFORM) ;
00087       M.ApplyOp(M.Ext.BC[i] , &M , i , INVERSEPOLYQUAD)  ;
00088     }
00089     
00090     // true solution
00091     M1.SetFunction(&FU,false,true) ;
00092     M1.Mul(0.5) ;
00093     
00094     // errors
00095     M.MMinus(&M1) ;
00096     l2er=sqrt( M.InnerProd(&M)/(e[0]*e[1]) ) ;
00097     
00098     // H1 error
00099     h1er=0; 
00100     for (i=0; i<2; i++) {
00101       M1.ApplyOp(M.Ext.BC[i] , &M , i , FD12) ;
00102       h1er += sqrt( M1.InnerProd(&M1)/(e[0]*e[1]) ) ;
00103     }
00104     
00105     printf("Error %d %e %e | %e %e \n",LL,l2er , l2era/l2er , h1er , h1era/h1er) ;
00106     l2era=l2er ;
00107     h1era=h1er ;
00108     
00109   }
00110 }
Generated at Mon Aug 19 10:02:32 2002 for AWFD by