| View previous topic :: View next topic |
| Author |
Message |
Axel Vogt
Guest
|
 Posted: Thu Nov 13, 2008 2:47 am Post subject: Question on Gauss quadrature with measure = cos |
 |
|
For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
(but did not find some searching for it), so used Maple to get it
from classical formulae for nodes and weights (using high precision
and converting to rationals appropriate for actual life).
In theory a rule of degree N should be exact for polynomials f up
to degree 2*N-1.
For N=4 that is almost ok, however beyond it becomes dirty and for
N=8 it begins to loose 4 - 8 digits (10 places would be ok for me).
Any advice (beyond "do not do it"), what is the reason?
And suggestion for the integral (only some evaluations, 10 places
of exactness to be expected) are welcome as well, of course ... |
|
| |
|
 Back to top |
user923005
Guest
|
| Posted: Thu Nov 13, 2008 2:47 am Post subject: Re: Question on Gauss quadrature with measure = cos |
|
|
On Nov 12, 12:47 pm, Axel Vogt <&nore...@axelvogt.de> wrote:
[quote]For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
(but did not find some searching for it), so used Maple to get it
from classical formulae for nodes and weights (using high precision
and converting to rationals appropriate for actual life).
In theory a rule of degree N should be exact for polynomials f up
to degree 2*N-1.
For N=4 that is almost ok, however beyond it becomes dirty and for
N=8 it begins to loose 4 - 8 digits (10 places would be ok for me).
Any advice (beyond "do not do it"), what is the reason?
And suggestion for the integral (only some evaluations, 10 places
of exactness to be expected) are welcome as well, of course ...
[/quote]
First of all, let me recommend double exponential quadrature (as
opposed to Gaussian quadrature) by reason of its admitted excellence.
Second, for an oscillating function, there are special versions.
See:
http://www.kurims.kyoto-u.ac.jp/~ooura/intde.html
This is also quite interesting:
http://crd.lbl.gov/~dhbailey/mpdist/
Arprec in the above distribution has a *really* well written arbitrary
precision double exponential quadrature routine.
This web search:
http://www.google.com/search?rls=ig&hl=en&q=double+exponential+quadrature&aq=f&oqmay also prove revealing. |
|
| |
|
| Back to top |
user923005
Guest
|
| Posted: Thu Nov 13, 2008 2:47 am Post subject: Re: Question on Gauss quadrature with measure = cos |
|
|
On Nov 12, 1:30 pm, Axel Vogt <&nore...@axelvogt.de> wrote:
[quote]user923005 wrote:
On Nov 12, 12:47 pm, Axel Vogt <&nore...@axelvogt.de> wrote:
For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
(but did not find some searching for it), so used Maple to get it
from classical formulae for nodes and weights (using high precision
and converting to rationals appropriate for actual life).
In theory a rule of degree N should be exact for polynomials f up
to degree 2*N-1.
For N=4 that is almost ok, however beyond it becomes dirty and for
N=8 it begins to loose 4 - 8 digits (10 places would be ok for me).
Any advice (beyond "do not do it"), what is the reason?
And suggestion for the integral (only some evaluations, 10 places
of exactness to be expected) are welcome as well, of course ...
First of all, let me recommend double exponential quadrature (as
opposed to Gaussian quadrature) by reason of its admitted excellence.
Second, for an oscillating function, there are special versions.
See:
http://www.kurims.kyoto-u.ac.jp/~ooura/intde.html
This is also quite interesting:
http://crd.lbl.gov/~dhbailey/mpdist/
Arprec in the above distribution has a *really* well written arbitrary
precision double exponential quadrature routine.
This web search:
http://www.google.com/search?rls=ig&hl=en&q=double+exponential+quadra...
may also prove revealing.
Thank you, I am "aware" of that. May be my question was to vague
(and finally I want to fall back to usual double precision as my
working environment): it is "what is or may be the trap" for the
sketched approach (besides the desire of fast evaluation)?
[/quote]
Now, obviously you are not integrating polynomials, which is fine --
of course -- if you have a smooth continuous function and you do not
have an h (width) that is too large.
Since you are integrating over a domain of pi in with, you must
subdivide it to get an accurate answer. How wide was your interval
for integration after subdivision?
I generally aim for around 0.25 or less with a 52 knott Gaussian rule
(and smaller if the derivatives are large -- in this case we do not
know unless you know a lot about f()).
Here is a simplistic automatic integration routine that calculates cos
(x) * f(x) for f(x) = x*x:
#include <float.h>
#include <stdlib.h>
#include <math.h>
typedef long double ( * ldIntegralFunction_t)(long double, const long
double *);
class Integral
{
//////////////////////////////////////////////////////////////////////////
// PUBLIC CLASS INFORMATION
//////////////////////////////////////////////////////////////////////////
public:
//////////////////////////////////////////////////////////////////////////
// Construction
//////////////////////////////////////////////////////////////////////////
Integral( // Constructor requires no arguments.
); // This version performs no evaluations on construction
Integral( // Constructor requires three arguments.
long double ldX0, // Beginning of the interval to integrate.
long double ldX1, // Ending of the interval to integrate.
ldIntegralFunction_t ldF // The derivative function to
integrate.
);
Integral( // Constructor requires three arguments.
long double ldX0, // Beginning of the interval to integrate.
long double ldX1, // Ending of the interval to integrate.
ldIntegralFunction_t ldF, // The derivative function to
integrate.
long double *pldConst // Optional const array is supplied for
this one.
);
//////////////////////////////////////////////////////////////////////////
// Destruction
//////////////////////////////////////////////////////////////////////////
~Integral(); // Destructor
//////////////////////////////////////////////////////////////////////////
// Inquiry Methods
//////////////////////////////////////////////////////////////////////////
// Display the result of the integration
void vDisplayResultsToConsole(); // Display the results to the
console
long double ldGetIntegralResult(); // Get the results of the
integration
long double ldGetX0(); // Get the start of the interval of
integration
long double ldGetX1(); // Get the end of the interval of
integration
//////////////////////////////////////////////////////////////////////////
// Modifier Methods
//////////////////////////////////////////////////////////////////////////
// Set the Derivative Function
void vSetFprime( ldIntegralFunction_t ldF );
void vSetX0( long double ldX0 ); // Set the start of the interval
void vSetX1( long double ldX1 ); // Set the end of the interval
void vSetConstArr( long double *pldConst ); // Set the optional
constant array
//////////////////////////////////////////////////////////////////////////
// Special Modifier Method
//////////////////////////////////////////////////////////////////////////
void Evaluate(); // Evaluate the integral
//////////////////////////////////////////////////////////////////////////
// PRIVATE CLASS INFORMATION
//////////////////////////////////////////////////////////////////////////
private:
//////////////////////////////////////////////////////////////////////////
// Data Attributes
//////////////////////////////////////////////////////////////////////////
ldIntegralFunction_t ldFprime; // Derivative Function.
long double ldIntegralResult; // Result of the integration.
long double ldX0; // Beginning of the interval to integrate.
long double ldX1; // Ending of the interval to integrate.
long double *pldConstArray; // Optional array of constants used by
function.
};
#include <float.h>
// internal prototype:
static long double qxrul(
ldIntegralFunction_t ldF,
long double ldA,
long double ldB,
int iRuleChoice,
long double *pldConstArray,
long double *pldVectorOne,
long double *pldVectorTwo,
int *pI1,
int *pI2
);
#define HALF_KNOT_COUNT 21
#define ABS( a ) ( a ) > 0 ? ( a ) : -( a )
///////////////////////////////////////////////////////////////////////////
// Function ldIntegrate() performs numerical integration over a
supplied
// interval (input parameters ldX0 to ldX1) for the supplied function.
//
// This formula used and the associated constants are from:
//
// the "Encyclopedic Dictionary of Mathematics, Volume 1"
// "1 Abel to 304 Numerical Solution of Partial Differential
Equations"
// By Kiyosi Ito.
//
// Imprint Cambridge, MA : MIT Press, 1993.
//
// ISBN 0262590204 ( set : PB )
// 0262090260 ( set : HB )
//
// LCCN 93139017
// Article 299 ( XV.7 ) Numerical Integration, pp 1120-1.
//
///////////////////////////////////////////////////////////////////////////
// Constructor:
//////////////////////////////////////////////////////////////////////////
Integral::Integral()
{
// Set up the integration interval values;
ldX0 = 0.0e0L;
ldX1 = 0.0e0L;
ldFprime = ( ldIntegralFunction_t ) NULL;
pldConstArray = ( long double * ) NULL;
}
Integral::Integral(
long double ldStart,
long double ldEnd,
ldIntegralFunction_t ldF
)
{
// Set up the integration interval values;
ldX0 = ldStart;
ldX1 = ldEnd;
ldFprime = ldF;
pldConstArray = ( long double * ) NULL;
Evaluate();
}
Integral::Integral(
long double ldStart,
long double ldEnd,
ldIntegralFunction_t ldF,
long double *pldConstArrayIn
)
{
// Set up the integration interval values;
ldX0 = ldStart;
ldX1 = ldEnd;
ldFprime = ldF;
pldConstArray = pldConstArrayIn;
Evaluate();
}
void Integral::Evaluate()
{
// Empirically derived for long double stepsize:
const long double ldSmallEnoughInterval = 0.25e0L;
// If interval is too large, divide the integral into two equal
sections.
if ( ( ldX1 - ldX0 ) > ldSmallEnoughInterval )
{
long double ldMidpoint = ldX0 + ( ldX1 - ldX0 ) * 0.5e0L;
Integral i0( ldX0, ldMidpoint, ldFprime, pldConstArray );
Integral I1( ldMidpoint, ldX1, ldFprime, pldConstArray );
ldIntegralResult = i0.ldGetIntegralResult() +
I1.ldGetIntegralResult();
}
else
{
const long double ldTol = LDBL_EPSILON*2.0e0L;
int I1 = -1;
int I2 = -1;
long double pldVectorOne[HALF_KNOT_COUNT] = { 0 };
long double pldVectorTwo[HALF_KNOT_COUNT] = { 0 };
int iRuleChoice = 0;
long double ldResult0;
long double ldResult1;
ldResult0 = qxrul( ldFprime, ldX0, ldX1, iRuleChoice++,
pldConstArray, pldVectorOne, pldVectorTwo, &I1, &I2 );
ldResult1 = qxrul( ldFprime, ldX0, ldX1, iRuleChoice++,
pldConstArray, pldVectorOne, pldVectorTwo, &I1, &I2 );
if ( ABS( ( ldResult0 - ldResult1 )/ldResult1 ) > ldTol )
{
ldResult0 = ldResult1;
ldResult1 = qxrul( ldFprime, ldX0, ldX1, iRuleChoice++,
pldConstArray, pldVectorOne, pldVectorTwo, &I1, &I2 );
if ( ABS( ( ldResult0 - ldResult1 )/ldResult1 ) > ldTol )
{
ldResult1 = qxrul( ldFprime, ldX0, ldX1, iRuleChoice,
pldConstArray, pldVectorOne, pldVectorTwo, &I1, &I2 );
}
}
ldIntegralResult = ldResult1;
}
}
///////////////////////////////////////////////////////////////////////////
// Destructor:
//////////////////////////////////////////////////////////////////////////
Integral::~Integral()
{
}
//////////////////////////////////////////////////////////////////////////
// Inquiry Methods:
//////////////////////////////////////////////////////////////////////////
void Integral::vDisplayResultsToConsole() // Display the results to
the console
{
}
long double Integral::ldGetIntegralResult(void)
{
return ldIntegralResult;
}
long double Integral::ldGetX0() // Get the start of the interval of
integration
{
return ldX0;
}
long double Integral::ldGetX1() // Get the end of the interval of
integration
{
return ldX1;
}
//////////////////////////////////////////////////////////////////////////
// Modifier Methods:
//////////////////////////////////////////////////////////////////////////
// Set the Derivative Function
void Integral::vSetFprime( ldIntegralFunction_t ldF )
{
ldFprime = ldF;
}
void Integral::vSetX0( long double ldX ) // Set the start of the
interval
{
ldX0 = ldX;
}
void Integral::vSetX1( long double ldX ) // Set the end of the
interval
{
ldX1 = ldX;
}
void Integral::vSetConstArr( long double *pldConst ) // Set the end of
the interval
{
pldConstArray = pldConst;
}
/*
** These routines are derived from function QXRUL from:
** ALGORITHM 691, COLLECTED ALGORITHMS FROM ACM.
** THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
** VOL. 17, NO. 2, PP. HALF_KNOT_COUNT8-232. JUNE, 1991.
**
** This was a component of the SLATEC FORTRAN library created by:
**
** LAWRENCE LIVERMORE NATIONAL LABORATORY
**
** To compute I = integral of f() over ( ldA, ldB ), with error
estimate
**
** Parameters on entry:
** f() - long double function to be integrated
**
** ldA - long double lower limit of integration
**
** ldB - long double upper limit of integration
**
**
** iRuleChoice - int
** Selector for choice of local integration rule. Rule used with:
** 13 points if iRuleChoice = 0,
** 19 points if iRuleChoice = 1,
** 27 points if iRuleChoice = 2,
** 42 points if iRuleChoice = 3
**
**
** pldVectorOne - long double
** vector containing pI1 saved functional values of positive
abscissas.
**
** pldVectorTwo - long double
** vector containing pI2 saved functional values of negative
abscissas.
**
** pI1, pI2 - int number of elements in pldVectorOne, pldVectorTwo
respectively.
**
** On return
** ldY - long double approximation to the integral I.
**
**
** pldVectorOne - long double
** vector containing pI1 saved functional values of positive
abscissas.
**
** pldVectorTwo - long double
** vector containing pI2 saved functional values of negative
abscissas.
**
** pI1, pI2 - int number of elements in pldVectorOne, pldVectorTwo
respectively.
*/
static long double qxrul(
ldIntegralFunction_t ldF,
long double ldA,
long double ldB,
int iRuleChoice,
long double *pldConstArray,
long double *pldVectorOne,
long double *pldVectorTwo,
int *pI1,
int *pI2
)
{
/* Starting positions into the weights array: */
int iStart[] = { 0, 7, 17, 31 };
/* Lengths of the weight sequence to process: */
int iLen[] = { 7, 10, 14, HALF_KNOT_COUNT };
/* Symmetrical offsets +/- for unit interval */
const long double ldAbscissae[HALF_KNOT_COUNT] {
0.0000000e0L, 0.2500e0L, 0.5000e0L, 0.750000e0L, 0.87500e0L,
0.9375000e0L, 1.0000e0L, 0.3750e0L, 0.625000e0L, 0.96875e0L,
0.1250000e0L, 0.6875e0L, 0.8125e0L, 0.984375e0L, 0.18750e0L,
0.3125000e0L, 0.4375e0L, 0.5625e0L, 0.843750e0L, 0.90625e0L,
0.9921875e0L
};
const long double ldWeights[52] {
0.1303262173284849021810473057638590518409112513421e0L,
0.2390632866847646220320329836544615917290026806242e0L,
0.2630626354774670227333506083741355715758124943143e0L,
0.2186819313830574175167853094864355208948886875898e0L,
0.027578976466428368658596011976074715743366742067e0L,
0.105575010053845844336503487908666979130555049383e0L,
0.01571194260595182254168429283636656908546309467968e0L,
0.1298751627936015783241173611320651866834051160074e0L,
0.2249996826462523640447834514709508786970828213187e0L,
0.1680415725925575286319046726692683040162290325505e0L,
0.1415567675701225879892811622832845252125600939627e0L,
0.1006482260551160175038684459742336605269707889822e0L,
0.02510604860724282479058338820428989444699235030871e0L,
0.009402964360009747110031098328922608224934320397592e0L,
0.0554269923329587516840678369514364633827480535978e0L,
0.09986735247403367525720377847755415293097913496236e0L,
0.04507523056810492466415880450799432587809828791196e0L,
0.06300942249647773931746170540321811473310938661469e0L,
0.1261383225537664703012999637242003647020326905948e0L,
0.127386443358102827287870998185030736345352311788e0L,
0.08576500414311820514214087864326799153427368592787e0L,
0.07102884842310253397447305465997026228407227220665e0L,
0.05026383572857942403759829860675892897279675661654e0L,
0.004683670010609093810432609684738393586390722052124e0L,
0.1235837891364555000245004813294817451524633100256e0L,
0.1148933497158144016800199601785309838604146040215e0L,
0.0125257577422612263339147770259358530725452719807e0L,
0.123957239623183424219418967424381861904228081664e0L,
0.02501306413750310579525950767549691151739047969345e0L,
0.04915957918146130094258849161350510503556792927578e0L,
0.02259167374956474713302030584548274729936249753832e0L,
0.06362762978782724559269342300509058175967124446839e0L,
0.09950065827346794643193261975720606296171462239514e0L,
0.07048220002718565366098742295389607994441704889441e0L,
0.06512297339398335645872697307762912795346716454337e0L,
0.03998229150313659724790527138690215186863915308702e0L,
0.03456512257080287509832054272964315588028252136044e0L,
0.002212167975884114432760321569298651047876071264944e0L,
0.08140326425945938045967829319725797511040878579808e0L,
0.06583213447600552906273539578430361199084485578379e0L,
0.02592913726450792546064232192976262988065252032902e0L,
0.1187141856692283347609436153545356484256869129472e0L,
0.05999947605385971985589674757013565610751028128731e0L,
0.0550093798019804173691025798834610183906258148982e0L,
0.005264422421764655969760271538981443718440340270116e0L,
0.01533126874056586959338368742803997744815413565014e0L,
0.03527159369750123100455704702965541866345781113903e0L,
0.05000556431653955124212795201196389006184693561679e0L,
0.05744164831179720106340717579281831675999717767532e0L,
0.01598823797283813438301248206397233634639162043386e0L,
0.02635660410220884993472478832884065450876913559421e0L,
0.01196003937945541091670106760660561117114584656319e0L
};
long double ldY = 0.0e0L;
long double ldOffset;
long double ldMidpoint;
long double ldHalfInterval;
int iKnots;
int i;
iKnots = iLen[iRuleChoice];
ldHalfInterval = ( ldB - ldA ) * 0.5;
ldMidpoint = ldA + ldHalfInterval;
for ( i = 0; i < iKnots; i++ )
{
if ( i >= *pI1 || i >= *pI2 )
{
ldOffset = ldHalfInterval * ldAbscissae[i];
if ( i >= *pI1 )
{
pldVectorOne[i] = ( *ldF )( ldMidpoint + ldOffset,
pldConstArray );
}
if ( i >= *pI2 )
{
pldVectorTwo[i] = ( *ldF )( ldMidpoint - ldOffset,
pldConstArray );
}
}
ldY += ( pldVectorOne[i] + pldVectorTwo[i] ) * ldWeights[iStart
[iRuleChoice] + i];
}
ldY *= ldHalfInterval;
if ( *pI1 <= iKnots )
{
*pI1 = iKnots;
}
if ( *pI2 <= iKnots )
{
*pI2 = iKnots;
}
return ldY;
}
static const long double ldPiOver2 1.57079632679489661923132169163975144209858469968755291e0L;
long double l( long double ldX, long double *pldNull )
{
return cosl(ldX) * ldX * ldX;
}
#include <stdio.h>
int main( int argc, char **argv )
{
Integral CosXLogX( -ldPiOver2, ldPiOver2, ( ldIntegralFunction_t )
l );
printf( "Integrated value = %.20Lf\n", CosXLogX.ldGetIntegralResult
() );
/*
I := 0.5 * pi * pi * sin(0.5 * pi)
- 4.0 * sin(0.5 * pi)
+ 2.0 * pi cos(0.5 * pi)
I := 1/2 *pi *pi - 4
I ~ 0.93480220054467931
Calculated on a machine with 8 byte double precision:
0.93480220054467944
*/
return 0;
} |
|
| |
|
| Back to top |
Axel Vogt
Guest
|
| Posted: Thu Nov 13, 2008 3:30 am Post subject: Re: Question on Gauss quadrature with measure = cos |
|
|
user923005 wrote:
[quote]On Nov 12, 12:47 pm, Axel Vogt <&nore...@axelvogt.de> wrote:
For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
(but did not find some searching for it), so used Maple to get it
from classical formulae for nodes and weights (using high precision
and converting to rationals appropriate for actual life).
In theory a rule of degree N should be exact for polynomials f up
to degree 2*N-1.
For N=4 that is almost ok, however beyond it becomes dirty and for
N=8 it begins to loose 4 - 8 digits (10 places would be ok for me).
Any advice (beyond "do not do it"), what is the reason?
And suggestion for the integral (only some evaluations, 10 places
of exactness to be expected) are welcome as well, of course ...
First of all, let me recommend double exponential quadrature (as
opposed to Gaussian quadrature) by reason of its admitted excellence.
Second, for an oscillating function, there are special versions.
See:
http://www.kurims.kyoto-u.ac.jp/~ooura/intde.html
This is also quite interesting:
http://crd.lbl.gov/~dhbailey/mpdist/
Arprec in the above distribution has a *really* well written arbitrary
precision double exponential quadrature routine.
This web search:
http://www.google.com/search?rls=ig&hl=en&q=double+exponential+quadrature&aq=f&oq=
may also prove revealing.
[/quote]
Thank you, I am "aware" of that. May be my question was to vague
(and finally I want to fall back to usual double precision as my
working environment): it is "what is or may be the trap" for the
sketched approach (besides the desire of fast evaluation)? |
|
| |
|
| Back to top |
Peter Spellucci
Guest
|
| Posted: Thu Nov 13, 2008 12:59 pm Post subject: Re: Question on Gauss quadrature with measure = cos |
|
|
In article <6o0touF156dsU1@mid.individual.net>,
Axel Vogt <&noreply@axelvogt.de> writes:
[quote]For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
[/quote]
^^^ really a multiply here , no singularity involved
why not using the classical Gauss-Legendre for g(x)=f(x)*cos(x)?
tables of weights and nodes for this for example inside
quadpack.
[quote](but did not find some searching for it), so used Maple to get it
from classical formulae for nodes and weights (using high precision
and converting to rationals appropriate for actual life).
In theory a rule of degree N should be exact for polynomials f up
to degree 2*N-1.
For N=4 that is almost ok, however beyond it becomes dirty and for
N=8 it begins to loose 4 - 8 digits (10 places would be ok for me).
Any advice (beyond "do not do it"), what is the reason?
[/quote]
the classical formulae are known to be highly unstable under roundoff
maybe your precision wasn>t high enough? cos(x) on [-pi/2,pi/2]
is a legal weight function , hence it should work.
did you try Gautschi>s orthpol from netlib/toms?
[quote]
And suggestion for the integral (only some evaluations, 10 places
of exactness to be expected) are welcome as well, of course ...
[/quote]
depends much on what f(x) is like ...
hth
peter |
|
| |
|
| Back to top |
Greg Heath
Guest
|
| Posted: Thu Nov 13, 2008 4:59 pm Post subject: Re: Question on Gauss quadrature with measure = cos |
|
|
On Nov 12, 3:47 pm, Axel Vogt <&nore...@axelvogt.de> wrote:
[quote]For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
(but did not find some searching for it), so used Maple to get it
from classical formulae for nodes and weights (using high precision
and converting to rationals appropriate for actual life).
In theory a rule of degree N should be exact for polynomials f up
to degree 2*N-1.
For N=4 that is almost ok, however beyond it becomes dirty and for
N=8 it begins to loose 4 - 8 digits (10 places would be ok for me).
Any advice (beyond "do not do it"), what is the reason?
And suggestion for the integral (only some evaluations, 10 places
of exactness to be expected) are welcome as well, of course ...
[/quote]
Have you tried Filon>s formula?
Hope this helps.
Greg |
|
| |
|
| Back to top |
Axel Vogt
Guest
|
| Posted: Fri Nov 14, 2008 2:22 am Post subject: Re: Question on Gauss quadrature with measure = cos |
|
|
user923005 wrote:
[quote]On Nov 12, 1:30 pm, Axel Vogt <&nore...@axelvogt.de> wrote:
user923005 wrote:
On Nov 12, 12:47 pm, Axel Vogt <&nore...@axelvogt.de> wrote:
For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
(but did not find some searching for it), so used Maple to get it
from classical formulae for nodes and weights (using high precision
and converting to rationals appropriate for actual life).
In theory a rule of degree N should be exact for polynomials f up
to degree 2*N-1.
For N=4 that is almost ok, however beyond it becomes dirty and for
N=8 it begins to loose 4 - 8 digits (10 places would be ok for me).
Any advice (beyond "do not do it"), what is the reason?
And suggestion for the integral (only some evaluations, 10 places
of exactness to be expected) are welcome as well, of course ...
First of all, let me recommend double exponential quadrature (as
opposed to Gaussian quadrature) by reason of its admitted excellence.
Second, for an oscillating function, there are special versions.
See:
http://www.kurims.kyoto-u.ac.jp/~ooura/intde.html
This is also quite interesting:
http://crd.lbl.gov/~dhbailey/mpdist/
Arprec in the above distribution has a *really* well written arbitrary
precision double exponential quadrature routine.
This web search:
http://www.google.com/search?rls=ig&hl=en&q=double+exponential+quadra...
may also prove revealing.
Thank you, I am "aware" of that. May be my question was to vague
(and finally I want to fall back to usual double precision as my
working environment): it is "what is or may be the trap" for the
sketched approach (besides the desire of fast evaluation)?
Now, obviously you are not integrating polynomials, which is fine --
of course -- if you have a smooth continuous function and you do not
have an h (width) that is too large.
Since you are integrating over a domain of pi in with, you must
subdivide it to get an accurate answer. How wide was your interval
for integration after subdivision?
I generally aim for around 0.25 or less with a 52 knott Gaussian rule
(and smaller if the derivatives are large -- in this case we do not
know unless you know a lot about f()).
Here is a simplistic automatic integration routine that calculates cos
(x) * f(x) for f(x) = x*x:
[/quote]
.... detailed C++ code snipped from the original answer ...
#include <stdio.h>
[quote]int main( int argc, char **argv )
{
Integral CosXLogX( -ldPiOver2, ldPiOver2, ( ldIntegralFunction_t )
l );
printf( "Integrated value = %.20Lf\n", CosXLogX.ldGetIntegralResult
() );
/*
I := 0.5 * pi * pi * sin(0.5 * pi)
- 4.0 * sin(0.5 * pi)
+ 2.0 * pi cos(0.5 * pi)
I := 1/2 *pi *pi - 4
I ~ 0.93480220054467931
Calculated on a machine with 8 byte double precision:
0.93480220054467944
*/
return 0;
}
[/quote]
Thank you! I will have to try it (currently working in Maple and
would have to compile it into some DLL for using it from there).
Working with 14 Digits (for technical reasons) I get the following
zeros and Gauss weights for measure=cos (degree 4 resp. 8 and sorry
for not caring for short notations using symmetry):
X4 := [-809566/679921, -5651323/12864749,
5651323/12864749, 809566/679921];
W4 := [747658/3336707, 5245870/6760757,
5245870/6760757, 747658/3336707];
X8 := [-4183527/2902223, -2436169/2111567,
-2820555/3793778, -5034555/19628162,
5034555/19628162, 2820555/3793778,
2436169/2111567, 4183527/2902223];
W8 := [951833/34567318, 542319/3760774,
608902/1810783, 801071/1628207,
801071/1628207, 608902/1810783,
542319/3760774, 951833/34567318];
The the 4-pt rule Sum( W4[i]*f(X4[i]), i=1..4) calculates the
integral for f(x)=x^2 as .93480220054415 with following errors:
`error [abs,rel]` = 0.53e-12, 0.57e-12 for the 4-pt rule and
`error [abs,rel]` = 0.65e-12, 0.695e-12 for the 8-pt rule.
A typical function I have is exp(-(x+10))*sin( Pi/8*(x - 1/10))
(slowly exponential decay and a possible, mildly oscillating
factor) in the application I have in mind. For that I get
`error [abs,rel]` = -0.26e-10, 0.11e-5 with 4 points
`error [abs,rel]` = -0.15e-16, 0.64e-12 with 8 points
While the relative error is not that nice (due to the method
I guess) it gives acceptable results with only 4 evaluations
of the function f (which is quite expensive in my real case
and that is the reason why I want to avoid many steps and do
accept errors - even if not reasonable estimated).
So I seem to loose 2 - 4 decimal places that way (and have
not checked, what your much more precise solution would need). |
|
| |
|
| Back to top |
Axel Vogt
Guest
|
| Posted: Fri Nov 14, 2008 2:23 am Post subject: Re: Question on Gauss quadrature with measure = cos |
|
|
Peter Spellucci wrote:
[quote]In article <6o0touF156dsU1@mid.individual.net>,
Axel Vogt <&noreply@axelvogt.de> writes:
For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
^^^ really a multiply here , no singularity involved
why not using the classical Gauss-Legendre for g(x)=f(x)*cos(x)?
tables of weights and nodes for this for example inside
quadpack.
(but did not find some searching for it), so used Maple to get it
from classical formulae for nodes and weights (using high precision
and converting to rationals appropriate for actual life).
In theory a rule of degree N should be exact for polynomials f up
to degree 2*N-1.
For N=4 that is almost ok, however beyond it becomes dirty and for
N=8 it begins to loose 4 - 8 digits (10 places would be ok for me).
Any advice (beyond "do not do it"), what is the reason?
the classical formulae are known to be highly unstable under roundoff
maybe your precision wasn>t high enough? cos(x) on [-pi/2,pi/2]
is a legal weight function , hence it should work.
did you try Gautschi>s orthpol from netlib/toms?
And suggestion for the integral (only some evaluations, 10 places
of exactness to be expected) are welcome as well, of course ...
depends much on what f(x) is like ...
hth
peter
[/quote]
Thx, it helped to locate my problem: my 'exact' reference values
for testing had huge rational coefficients, so I have to care for
that (using different precisions for reference and result and only
after passing to a common precision that gives sense).
Yes, f(x)*cos(x) means multiplication. The classical method would
mean to 'write' that as polynomial, but having 0 in both endpoints
I think one would need higher order (and I want only few evaluation
points). So I was asking me, whether one could regard the cos as a
weight function. Comparing 4-point rules this seems to give better.
And it was interest in just experimenting and trying something else
(and was inspired by the Hurwitz-Zweifel method sketched in Kythe-
Schaferkotter, but I find that 'book' unenjoyable to read and do
not want to rely on that writings).
Typically my f is of exponential decay with a slowly oscillating
factor, something like f(x) = exp(-(x+10))*sin( Pi/8*(x - 1/10)). |
|
| |
|
| Back to top |
Axel Vogt
Guest
|
| Posted: Fri Nov 14, 2008 2:23 am Post subject: Re: Question on Gauss quadrature with measure = cos |
|
|
Greg Heath wrote:
[quote]On Nov 12, 3:47 pm, Axel Vogt <&nore...@axelvogt.de> wrote:
For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
(but did not find some searching for it), so used Maple to get it
from classical formulae for nodes and weights (using high precision
and converting to rationals appropriate for actual life).
In theory a rule of degree N should be exact for polynomials f up
to degree 2*N-1.
For N=4 that is almost ok, however beyond it becomes dirty and for
N=8 it begins to loose 4 - 8 digits (10 places would be ok for me).
Any advice (beyond "do not do it"), what is the reason?
And suggestion for the integral (only some evaluations, 10 places
of exactness to be expected) are welcome as well, of course ...
Have you tried Filon>s formula?
Hope this helps.
Greg
[/quote]
I think that is for f(x)*cos(w*x) for w << 1, but I have w=1. |
|
| |
|
| Back to top |
Greg Heath
Guest
|
| Posted: Fri Nov 14, 2008 1:52 pm Post subject: Re: Question on Gauss quadrature with measure = cos |
|
|
On Nov 13, 3:23 pm, Axel Vogt <&nore...@axelvogt.de> wrote:
[quote]Greg Heath wrote:
On Nov 12, 3:47 pm, Axel Vogt <&nore...@axelvogt.de> wrote:
For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
(but did not find some searching for it), so used Maple to get it
from classical formulae for nodes and weights (using high precision
and converting to rationals appropriate for actual life).
In theory a rule of degree N should be exact for polynomials f up
to degree 2*N-1.
For N=4 that is almost ok, however beyond it becomes dirty and for
N=8 it begins to loose 4 - 8 digits (10 places would be ok for me).
Any advice (beyond "do not do it"), what is the reason?
And suggestion for the integral (only some evaluations, 10 places
of exactness to be expected) are welcome as well, of course ...
Have you tried Filon>s formula?
Hope this helps.
Greg
I think that is for f(x)*cos(w*x) for w << 1, but I have w=1
[/quote]
No.
You are wrong.
http://groups.google.com/group/sci.math.num-analysis/msg/02c1c8667c8680f0
Hope this helps.
Greg |
|
| |
|
| Back to top |
Peter Spellucci
Guest
|
| Posted: Fri Nov 14, 2008 2:00 pm Post subject: Re: Question on Gauss quadrature with measure = cos |
|
|
In article <6o3gmuF1migqU2@mid.individual.net>,
Axel Vogt <&noreply@axelvogt.de> writes:
[quote]Peter Spellucci wrote:
In article <6o0touF156dsU1@mid.individual.net>,
Axel Vogt <&noreply@axelvogt.de> writes:
For Int( f(x)*cos(x), x = -Pi/2 ... Pi/2) i want a Gauss quadrature
[/quote]
[quote]Yes, f(x)*cos(x) means multiplication. The classical method would
mean to 'write' that as polynomial, but having 0 in both endpoints
[/quote]
??? if f is smooth, so f(.)*cos(.) is and
setting g(x)=f(x)*cos(x) would also not mean growing derivatives
of g relative to those of f, hence I cannot see why not
direct use of
sum w_i*g(t_i)
with the weights and nodes of gauss-legendre, linearly transformed to
your interval, (means multiplying both by pi/2) should have problems.
[quote]I think one would need higher order (and I want only few evaluation
points). So I was asking me, whether one could regard the cos as a
weight function. Comparing 4-point rules this seems to give better.
And it was interest in just experimenting and trying something else
(and was inspired by the Hurwitz-Zweifel method sketched in Kythe-
Schaferkotter, but I find that 'book' unenjoyable to read and do
not want to rely on that writings).
Typically my f is of exponential decay with a slowly oscillating
factor, something like f(x) = exp(-(x+10))*sin( Pi/8*(x - 1/10)).
I tried it myself:[/quote]
Formula: Gauss-Kronrod-Paar with 7 resp. 15 knots,
g(x)=exp(-(x+10))*sin((pi/8)*(x-0.1))*cos(x) on [-pi/2,pi/2]
# f-values = 15 (d.h. only one step,
no adaptivity needed)
Integral = -0.2343155286749262E-04
Errorestimate = 0.6889301389209115E-14
since the error estimate is simply the difference between
those two formulae, 7 knots are sufficient anyway
Maybe in general it is advantegeous to factor out the term exp(-10) or
similars
hth
peter |
|
| |
|
| Back to top |
|
