stration55125.cn

HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULESBRADLEYK.
ALPERTSIAMJ.
SCI.
COMPUT.
c1999SocietyforIndustrialandAppliedMathematicsVol.
20,No.
5,pp.
1551–1584Abstract.
Anewclassofquadraturerulesfortheintegrationofbothregularandsingularfunctionsisconstructedandanalyzed.
Foreachrulethequadratureweightsarepositiveandtheclassincludesrulesofarbitrarilyhigh-orderconvergence.
Thequadraturesresultfromalterationstothetrapezoidalrule,inwhichasmallnumberofnodesandweightsattheendsoftheintegrationintervalarereplaced.
Thenewnodesandweightsaredeterminedsothattheasymptoticexpansionoftheresultingrule,providedbyageneralizationoftheEuler–Maclaurinsummationformula,hasaprescribednumberofvanishingterms.
Thesuperiorperformanceoftherulesisdemonstratedwithnumericalexamplesandapplicationtoseveralproblemsisdiscussed.
Keywords.
Euler–Maclaurinformula,Gaussianquadrature,high-orderconvergence,numericalintegration,positiveweights,singularityAMSsubjectclassications.
41A55,41A60,65B15,65D32PII.
S10648275973251411.
Introduction.
Recentadvancesinalgorithmsforthenumericalsolutionofintegralequationshavestimulatedrenewedinterestinintegralequationformulationsofproblemsinpotentialtheory,wavepropagation,andotherapplicationareas.
Fastalgorithms,includingthosebyRokhlin[1],[2],GreengardandRokhlin[3],HackbuschandNowak[4],Beylkin,Coifman,andRokhlin[5],Alpert,Beylkin,Coifman,andRokhlin[6],andKelley[7],havegenerallyreducedthecomputationalcomplexitytoO(n)orO(nlogn)operations,withnunknowns,fortheapplicationofanintegraloperatororitsinverse.
Theappearanceofthesefastalgorithmshasincreasedtheurgencyofdevelopingaccuratequadraturesforthediscretizationofintegraloperators.
Suchquadraturesmusteectivelytreatthekernelsingularitiesoftheoperatorsandthevaryinglocationofthesingularitieswithaxedsetofdensitynodes,andmustallowtheapplicationofthefastalgorithms.
Inthispaper,wedevelopquadraturerules,basedonalterationstothetrapezoidalrule,thatobeytheseconstraints.
Thewell-knownEuler–Maclaurinsummationformulaprovidesanasymptoticex-pansionforthetrapezoidalruleappliedtoregularfunctions.
Whiletheconstanttermoftheexpansionisanintegral,theothertermsdependontheintegrand'sderivativesattheendpointsoftheintervalofintegration.
Thisexpansionisoftenusedto"cor-rect"thetrapezoidalruletoaquadraturewithhigh-orderconvergence,throughtheuseofeitherknownderivativevaluesortheirnite-dierenceapproximations.
Agen-eralizationoftheEuler–MaclaurinformulabyNavot[8],forintegrablefunctionswithsingularitiesoftheformxγ,canalsobeusedtocorrectthetrapezoidalruletoachievehigh-orderconvergence(providedtheexponentγisknown).
ReceivedbytheeditorsJuly29,1997;acceptedforpublication(inrevisedform)February20,1998;publishedelectronicallyApril27,1999.
ThisresearchwassupportedinpartbyDefenseAdvancedResearchProjectsAgencycontract9760400andinpartbytheUniversityofColoradoPrograminAppliedMathematics,Boulder,wheretheauthorwaswelcomedduringthegovernmentfurloughsof1995–1996.
ContributionofU.
S.
Government;notsubjecttocopyrightintheUnitedStates.
http://www.
siam.
org/journals/sisc/20-5/32514.
htmlNationalInstituteofStandardsandTechnology,325Broadway,Boulder,CO80303(alpert@boulder.
nist.
gov).
15511552BRADLEYK.
ALPERTWederivenewquadratures,basedontheEuler–Maclaurinformulaanditsgener-alization,ofarbitrarilyhigh-orderconvergence,forregularorsingularfunctionswithpowerorlogarithmsingularity.
Eachquadratureisconstructedbychangingthetrape-zoidalrule:afewofthenodesandweightsattheintervalendpointsarereplacedwithnewnodesandweightsdeterminedsoastoannihilateseveraltermsintheasymptoticexpansion.
Thenodesalwaysliewithintheintervalofintegrationandtheweightsarealwayspositive.
Foraregularfunctionf:[0,1]→R,weapproximate10f(x)dxwiththequadra-tureTn(f)=hw1f(x1h)+w2f(x2h)wjf(xjh)(1)+f(ah)+f(ah+h)f(1ah)+wjf(1xjh)w1f(1x1h).
Therearen"internal"nodeswithspacingh=1/(n+2a1)andj"endpoint"nodesateachend,withtheendpointnodesx1,xjandweightsw1,wjchosensothattheasymptoticexpansionofTnasn→∞has2jvanishingtermsandTn(f)=10f(x)dx+O(h2j+1)(2)(Theorem3.
1andCorollary3.
2).
Theparametersaandj,andthenodesx1,xjandweightsw1,wj,areindependentofn.
Thenodesandweightsaredeterminedby2jnonlinearequations,whichhaveauniquesolution,with00,i=1,j,(3)providedaissucientlylarge(Theorem4.
7).
Forintegrandsthataresingularatoneendpoint,Tnisalteredsothatthenodesandweightsatthatenddierfromthoseattheotherendanddependonthesingularity(Theorem3.
4andCorollary3.
6;Theorem3.
7andCorollary3.
8).
Forimproperintegralsinwhichtheintegrandisoscillatoryandslowlydecaying,TniscombinedwithGauss–Laguerrequadraturetogiveruleswithhigh-orderconvergence(Theorem3.
9andCorollary3.
10).
Severalauthorshavestudiedtheproblemstreatedhere.
Ithasbeenobservedthatendpointcorrectionscanbederivedforsingularintegrands;Rokhlin[9]implementedsuchaschemeforintegrandswithaknownsingularityatanintervalendpoint.
Hede-rivedcorrectionstothetrapezoidalrulebyplacingadditionalquadraturenodesneartheendpoint,withthecorrespondingweightsdeterminedsothatlow-orderpolyno-mialsandthesingularitytimeslow-orderpolynomialswereintegratedexactly.
Heshowedthatunderfairlygeneralconditions,theseweightshadlimitingvalues(uptoscale)asthenumberofnodesinthetrapezoidalruleincreasedwithoutboundandthattheselimitingweightscouldbeusedtoformquadratureruleswithgoodconver-gence.
Unfortunately,theorderofconvergenceoftheserulesisrestrictedinpracticebythefactthattheweightsincreaseinmagnituderapidlyastheorderincreases.
EortsbyStarr[10]andsubsequentlybyAlpert[11]reducedthegrowthinsizeoftheweightswithorder,primarilybyusingmoreweightsthanthenumberofequationssatisedandminimizingtheirsumofsquares.
Inanotherapproach,Kress[12]usesallquadratureweightsinthequadraturerule,ratherthanafewneartheendpoints,HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1553tohandlethesingularity.
Morerecently,KapurandRokhlin[13]successfullycon-structedrulesofarbitraryorderbyseparatingtheintegrand'sregularandsingularpartsandallowingsomequadraturenodestolieoutsidetheintervalofintegration.
Thepresentapproachdoesnotsuerfromlimitationsonorderofconvergence,separationoftheintegrandintoparts,orquadraturenodesoutsidetheintervalofconvergence.
Ontheotherhand,thequadraturenodesneartheintervalendpointsarenotequispaced.
Also,theequationsforthenodesx1,xjandweightsw1,wj,inadditiontobeingnonlinear,arepoorlyconditioned;theconditioningdeterioratesrapidlywithincreasingorder.
Nevertheless,weareabletouseanalgorithmdevelopedrecentlybyMa,Rokhlin,andWandzura[14]forcomputinggeneralizedGaussianquadraturestoobtainaccuratequadraturenodesandweights.
Theauthorwouldalsoliketocreditthatpaperforinspiringthepresentwork.
Thepaperisorganizedaroundsection3,wherethenewasymptoticexpansionsarederivedandthequadraturesdened,andsection4,whereitisshownthattheequationsdeningthequadraturesactuallyhavesolutions,whichareunique.
Thesesectionsareprecededbymathematicalpreliminariesandfollowedbyadiscussionofthecomputationofthequadraturenodesandweights.
Numericalexamplesarepresentedinsection6andweconcludewithsomeapplicationsandasummary.
2.
Mathematicalpreliminaries.
Thematerialinthissection,whichisfoundinstandardreferences,isusedinthesubsequentdevelopment.
2.
1.
Bernoullipolynomials.
TheBernoullipolynomialsaredenedbygener-atingfunction(see,forexample,[15,(23.
1.
1)])textet1=∞n=0Bn(x)tnn!
fromwhichB0(x)=1,B1(x)=x12,B2(x)=x2x+16,B3(x)=x332x2+12x.
TheBernoullipolynomialssatisfythedierenceformulaBn(x+k)Bn(x)n=k1i=0(x+i)n1,n=1,2,(4)thedierentiationformulaBn(x)=nBn1(x),n=1,2,(5)andtheexpansionformulaBn(x+h)=nr=0nrBr(x)hnr,n=0,1,(6)2.
2.
Euler–Maclaurinsummationformula.
Forafunctionf∈Cp(R),p≥1,theEuler–Maclaurinsummationformula(see,forexample,[15,(23.
1.
30)])canbederivedbyrepeatedintegrationbyparts.
Werstconsidertheinterval[c,c+h]and1554BRADLEYK.
ALPERTapply(5)toobtainc+hcf(x)dx=h10B0(1x)f(c+xh)dx(7)=hB1(1x)1!
f(c+xh)10+h210B1(1x)1!
f(c+xh)dx.
.
.
=p1r=0hr+1Br+1(1x)(r+1)!
f(r)(c+xh)10+hp+110Bp(1x)p!
f(p)(c+xh)dx=hf(c)+f(c+h)2p1r=1hr+1Br+1(r+1)!
f(r)(c+h)f(r)(c)+hp+110Bp(1x)p!
f(p)(c+xh)dx,wherewehaveusedB1(0)=B1(1)=12,Bn(0)=Bn(1)=Bn,n=1.
ToderivetheEuler–Maclaurinformulafortheinterval[a,b],weleth=(ba)/nandc=a+ihin(7),sumoveri=0,1,n1,andrearrangetermstoobtainhf(a)2+f(a+h)f(bh)+f(b)2(8)=baf(x)dx+p1r=1hr+1Br+1(r+1)!
f(r)(b)f(r)(a)hp+110Bp(1x)p!
n1i=0f(p)(a+ih+xh)dx.
Theexpressionontheleft-handsideof(8)isthewell-knowntrapezoidalrule.
Eval-uationoftheexpressionontheright-handsideof(8)issimpliedbythefactthatB2r+1=0forr≥1.
2.
3.
GeneralizedRiemannzeta-function.
ThegeneralizedRiemannzeta-functionisdenedbytheformulaζ(s,v)=∞n=01(v+n)s,Re(s)>1,v=0,1,.
.
.
.
Thisfunctionhasacontinuationthatisanalyticintheentirecomplexs-plane,withtheexceptionofs=1,whereithasasimplepole.
Inwhatfollows,weshallbeconcernedprimarilywithrealsandv,withs0.
WewillusetheHYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1555followingrepresentationderivedfromPlana'ssummationformula(see,forexample,[16,section1.
10(7)]):ζ(s,v)=v1ss1+vs2+2v1s∞0sin(sarctant)(1+t2)s/2dte2πvt1,Re(v)>0.
(9)Equation(9)canbeusedtoderivetheasymptoticexpansionofζasv→∞.
WetreattheintegralasasumofLaplaceintegrals,eachwithanasymptoticexpansiongivenbyWatson'slemma(see,forexample,[17,p.
263]),andobtainζ(s,v)=v1ss1+vs2+1s1pr=1s+2r22rB2rvs+2r1+O(vs2p1),(10)asv→∞,withs∈C,s=1,andpanarbitrarypositiveinteger.
Equation(10)isaslightgeneralizationof[16,section1.
18(9)].
ThereisadirectconnectionbetweentheBernoullipolynomialsandζ,Bn(x)n=ζ(1n,x),n=1,2,andgeneralizationsofthedierenceanddierentiationformulaehold:ζ(s,v)ζ(s,v+k)=k1i=01(v+i)s,(11)ζ(s,v)v=sζ(s+1,v).
(12)2.
4.
OrthogonalpolynomialsandGaussianquadrature.
Supposethatωisapositivecontinuousfunctionontheinterval(a,b)andωisintegrableon[a,b].
Wedenetheinnerproductwithrespecttoωofreal-valuedfunctionsfandgbytheintegral(f,g)=baf(x)g(x)ω(x)dx.
Thereexistpolynomialsp0,p1,ofdegree0,1,respectively,suchthat(pn,pm)=0forn=m(orthogonality);theyareuniqueuptothechoiceofleadingcoecients.
Withleadingcoecientsone,theycanbeobtainedrecursivelybytheformulae(see,forexample,[18,p.
143])p0(x)=1,(13)pn+1(x)=(xδn+1)pn(x)γn+12pn1(x),n≥0,(14)wherep1(x)=0andδn,γnaredenedbytheformulaeδn+1=(xpn,pn)/(pn,pn),n≥0,γn+12=0,n=0,(pn,pn)/(pn1,pn1),n≥1.
Thezerosxn1xnnofpnaredistinctandlieintheinterval(a,b).
Thereexistpositivenumbersωn1ωnnsuchthatbaf(x)ω(x)dx=ni=1ωnif(xni)(15)1556BRADLEYK.
ALPERTwheneverfisapolynomialofdegreelessthan2n.
TheseChristoelnumbersaregivenbytheformula(see,forexample,[19,p.
48])ωni=(pn1,pn1)pn1(xni)pn(xni),i=1,n.
(16)Moreover,ifω(x)=(bx)τ(x)withτintegrableon[a,b],then,withthedenitionxnn+1=b,thereexistpositivenumbersτn1τnn+1suchthatbaf(x)τ(x)dx=n+1i=1τnif(xni)(17)wheneverfisapolynomialofdegreelessthanorequalto2n.
ThesemodiedChristof-felnumbersaregivenbytheformulaτni=ωnibxni,i=1,n,bapn(x)pn(b)τ(x)dx,i=n+1,(18)whereωniisgivenby(16).
Thesummationin(15)isthen-nodeGaussianquadraturewithrespecttoω,whilethatin(17)isan(n+1)-nodeGauss–Radauquadraturewithrespecttoτ.
3.
HybridGauss-trapezoidalquadraturerules.
Inthissectionweintro-ducenewquadraturerulesforregularintegrands,singularintegrandswithapowerorlogarithmicsingularity,andimproperintegralsanddeterminetheirrateofconver-genceasthenumberofquadraturenodesincreases.
Fornotationalconveniencewegenerallyconsiderquadraturesoncanonicalinter-vals,primarily[0,1].
Itisunderstoodthatthesearereadilytransformedtoquadra-turesonanyniteinterval[a,b]bytheappropriatelineartransformationofthenodesandweights.
3.
1.
Regularintegrands.
Forj,npositiveintegersanda∈R+={x∈R|x>0},wedenealinearoperatorTjanonC([0,1]),dependingonnodesx1,xjandweightsw1,wj,bytheformulaTjan(f)=hji=1wif(xih)+hn1i=0f(ah+ih)+hji=1wif(1xih),(19)whereh=(n+2a1)1ischosensothatah+(n1)h=1ah.
Theorem3.
1.
Supposef∈Cp([0,1]).
TheasymptoticexpansionofTjan(f)asn→∞isgivenbytheformula(20)Tjan(f)=10f(x)dx+p1r=0hr+1f(r)(0)+(1)rf(r)(1)r!
ji=1wixirBr+1(a)r+1+O(hp+1).
HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1557Proof.
WeapplytheEuler–Maclaurinformula(8)ontheinterval[ah,1ah]toobtain(21)hn1i=0f(ah+ih)=hf(ah)+f(1ah)2+1ahahf(x)dx+p1r=1hr+1Br+1(r+1)!
f(r)(1ah)f(r)(ah)+O(hp+1).
Wenowcombine(19)and(21),theequality1ahahf(x)dx=10f(x)dxah0f(x)dx11ahf(x)dx,Taylorexpansionofallquantitiesabouth=0,theBernoullipolynomialexpansionformula(6),anddierenceformula(4)toobtain(20).
Corollary3.
2.
Supposethenodesx1,xjandweightsw1,wjsatisfytheequationsji=1wixir=Br+1(a)r+1,r=0,1,2j1.
(22)ThenTjanisaquadraturerulewithconvergenceoforder2j+1forf∈Cp([0,1])withp≥2j.
Moreover,Tjan(f)10f(x)dxh2j+1f(2j)(0)+f(2j)(1)(2j)!
ji=1wixi2jB2j+1(a)2j+1(23)asn→∞,providedf(2j)(0)+f(2j)(1)=0.
Corollary3.
3.
Supposexj=a1andtheremainingnodesx1,xj1andweightsw1,wjsatisfytheequationsji=1wixir=Br+1(a)r+1,r=0,1,2j2.
(24)ThenTjanisaquadraturerulewithconvergenceoforder2jforf∈Cp([0,1])withp≥2j1.
Moreover,Tjan(f)10f(x)dxh2jf(2j1)(0)f(2j1)(1)(2j1)!
ji=1wixi2j1B2j(a)2j(25)asn→∞,providedf(2j1)(0)f(2j1)(1)=0.
Weshallseebelowthat(22)hasasolutionwiththenodesandweightsallpositiveifaissucientlylargeandthatnumericalsolutionof(22)isequivalenttocomputingtherootsofaparticularpolynomial.
Thisstatementholdsfor(24)aswell.
3.
2.
Singularintegrands.
Forj,k,npositiveintegersanda,b∈R+,wedenealinearoperatorSjkabnonC((0,1]),dependingonnodesv1,vj,x1,xkandweightsu1,uj,w1,wk,bytheformulaSjkabn(g)=hji=1uig(vih)+hn1i=0g(ah+ih)+hki=1wig(1xih),(26)1558BRADLEYK.
ALPERTwhereh=(n+a+b1)1ischosensothatah+(n1)h=1bh.
Thefollowingtheorem,whichfollowsfromageneralizationoftheEuler–MaclaurinformuladuetoNavot[8]orafurthergeneralizationduetoLyness[20],presentsasomewhatdierentproofthantheearlierones.
Theorem3.
4.
Supposeg(x)=xγf(x),whereγ>1andf∈Cp([0,1]).
TheasymptoticexpansionofSjkabn(g)asn→∞isgivenbytheformula(27)Sjkabn(g)=10g(x)dx+p1r=0hγ+r+1f(r)(0)r!
ji=1uiviγ+r+ζ(γr,a)+p1r=0hr+1(1)rg(r)(1)r!
ki=1wixirBr+1(b)r+1+O(hp+1+min{0,γ}).
Proof.
Forc∈R+,wedenepolynomialspc0,pc1,inanalogywiththeBernoullipolynomials,bytheformulapcn(x)=nr=0nrζ(γr,c)(1cx)nr.
(28)Dierentiating,weverifythatddxpcn(x)=npcn1(x),n=1,2,.
.
.
and,combiningtheζdierenceformula(11)with(28),weobtainpcn(1)pc+1n(0)=cγ,n=0,0,n≥1.
Additionally,wedenefunctionsqc0,qc1,.
.
.
bytheformulaqcn(x)=(x+c)γ,n=0,(1)n(x+c)γ+n(γ+1)γ+n)pc+1n1(x)(n1)!
,n=1,2,(29)andobservethatddxqcn(x)=qcn1(x),n=1,2,HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1559Withthesedenitions,theprooffollowsthatoftheEuler–Maclaurinformula:1bhahxγf(x)dx=n2i=0h1+γ10qa+i0(x)f(ah+ih+xh)dx(30)=n2i=0p1r=0hγ+r+1qa+ir+1(x)f(r)(ah+ih+xh)10+hγ+p+110qa+ip(x)f(p)(ah+ih+xh)dx=hn2i=1(ah+ih)γf(ah+ih)p1r=0hγ+r+1qa+n2r+1(1)f(r)(1bh)qar+1(0)f(r)(ah)+hγ+p+110n2i=0qa+ip(x)f(p)(ah+ih+xh)dx.
Taylorexpansionoff(r)(ah)abouth=0,thedenitions(28)and(29)forpcnandqcn,andthebinomialtheoremcombinetoyield(31)p1r=0hγ+r+1qar+1(0)f(r)(ah)=p1k=0hγ+k+1f(k)(0)k!
aγ+k+1γ+k+1+ζ(γk,a+1)+O(hγ+p+1).
Likewise,Taylorexpansionoff(r)(1bh)abouth=0,thedenitionsforpcnandqcn,theasymptoticexpansion(10)forζ,theBernoullipolynomialexpansionformula(6),andthebinomialtheoremcombinetoyield(32)p1r=0hγ+r+1qa+n2r+1(1)f(r)(1bh)=p1k=0hk+1(1)kg(k)(1)k!
bk+1k+1Bk+1(b+1)k+1+O(hp+1).
Wenowcombine(26)and(30)–(32),theequality1bhahg(x)dx=10g(x)dxah0g(x)dx11bhg(x)dx,expansionofthelattertwointegralsabouth=0,andthedierenceformula(4)fortheBernoullipolynomialsand(11)forζtoobtain(27).
Corollary3.
5.
Supposethenodesu1,ujandweightsv1,vjsatisfytheequationsji=1uiviγ+r=ζ(γr,a),r=0,1,2j1,(33)1560BRADLEYK.
ALPERTandthenodesx1,xjandweightsw1,wjsatisfytheequationsji=1wixir=Br+1(b)r+1,r=0,1,2j1.
(34)ThenSjjabnisaquadraturerulewithconvergenceoforder2j+1+min{0,γ}forg,whereg(x)=xγf(x),withf∈Cp([0,1])andp≥2j.
Moreover,(35)Sjjabn(g)10g(x)dxhγ+2j+1f(2j)(0)(2j)!
ji=1uiviγ+2j+ζ(γ2j,a),γ0,g(2j)(1)=0,asn→∞.
Corollary3.
6.
Supposethenodesu1,ujandweightsv1,vjsatisfytheequationsji=1uiviγ+r=ζ(γr,a),r=0,1,j1,ji=1uivir=Br+1(a)r+1,r=0,1,j1,(36)andthenodesx1,xkandweightsw1,wksatisfytheequationski=1wixir=Br+1(b)r+1,r=0,1,2k1.
(37)ThenSjkabnisaquadraturerulewithconvergenceofordermin{j+1,γ+j+1,2k+1}forg,whereg(x)=xγφ(x)+ψ(x),withφ,ψ∈Cp([0,1])andp≥min{j,2k}.
InCorollaries3.
5and3.
6,anevennumberofconstraintsonthenodesandweightsatbothendsoftheintervalareconsidered.
Clearly,thereareanalogousquadraturerulesarisingfromanoddnumberofconstraintsatoneorbothends;thesearesimilar,andexplicitpresentationofthemisomitted.
Wenowconsideradierenttypeofsingularity.
Theorem3.
7.
Supposeg(x)=f(x)logx,wheref∈Cp([0,1]).
TheasymptoticexpansionofSjkabn(g)asn→∞isgivenbytheformulaSjkabn(g)=10g(x)dx(38)+p1r=0hr+1f(r)(0)r!
ji=1uivirlog(vjh)ζ(r,a)Br+1(a)r+1logh+p1r=0hr+1(1)rg(r)(1)r!
ki=1wixirBr+1(b)r+1+O(hp+1logh),HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1561whereζdenotesthederivativeofζwithrespecttoitsrstargument.
Proof.
ThisasymptoticexpansionisderivedfromthatofTheorem3.
4bydier-entiating(27)withrespecttoγandevaluatingtheresultatγ=0.
Corollary3.
8.
Supposethenodesu1,ujandweightsv1,vjsatisfytheequationsji=1uivirlogvi=ζ(r,a),r=0,1,j1,ji=1uivir=Br+1(a)r+1,r=0,1,j1,(39)andthenodesx1,xkandweightsw1,wksatisfytheequationski=1wixir=Br+1(b)r+1,r=0,1,2k1.
(40)ThenSjkabnisaquadraturerulewitherroroforderO(min{hj+1logh,h2k+1})forg,whereg(x)=φ(x)logx+ψ(x),withφ,ψ∈Cp([0,1])andp≥min{j,2k}.
3.
3.
Improperintegrals.
Forj,npositiveintegers,wedenealinearoperatorRjnonC([n,∞)),dependingonnodesx1,xjandweightsw1,wj,bytheformulaRjn(g)=jk=1wkg(n+xk).
(41)Theorem3.
9.
Supposeg(x)=eiγxf(x),whereγ∈R,γ=0,andf∈Cp([1,∞)),andthatthereexistpositiveconstantsβ,α0,αp,suchthatf(r)(x)1,r=0,p.
(42)TheasymptoticexpansionofRjn(g)asn→∞isgivenbytheformula(43)Rjn(g)=∞ng(x)dx+eiγnp1r=0f(r)(n)r!
jk=1wkxkreiγxkr!
iγr+1+O(np).
Proof.
Weintegratebypartsrepeatedlytoobtain∞neiγxf(x)dx=eiγnp1r=0iγr+1f(r)(n)+∞neiγxiγpf(p)(x)dx,(44)andin(41)wecomputetheTaylorexpansionoffaboutxk=0togetRjn(g)=eiγnjk=1eiγxkwkp1r=0f(r)(n)r!
xkr+f(p)(n+ξik)p!
xkp,(45)1562BRADLEYK.
ALPERTwhereξkliesbetween0andxkfork=1,j.
Nowcombining(42),(44),and(45),weobtain(43).
Example.
Thefunctionfdenedbytheformulaf(x)=∞r=0arxβ+r,(46)with|ar|aforsomea∈R.
Supposefurtherthatv1,vjaretherootsoftheLaguerrepolynomialLjofdegreej,thatcoecientsu1,ujsatisfytheequationsjk=1ukvkrevk=r!
,r=0,1,j1,(47)andthattheoperatorRjnisdenedwithnodesxk=(i/γ)vkandweightswk=(i/γ)ukfork=1,j.
SupposenallythatTjamisdenedtobethequadratureruleTjamwithnodesandweightssatisfying(22)buttranslatedandscaledtotheinterval[1,n].
Thenforp≥2j,theexpressionTja(n1)n(g)+Rjn(g)isanapproximationfortheintegral∞1g(x)dx,whereg(x)=eiγxf(x),witherroroforderO(n2j)asn→∞.
Proof.
ThisresultisjustacombinationofthequadratureruleofCorollary3.
2,fortheinterval[1,n],withtheasymptoticexpansionofTheorem3.
9,fortheinterval[n,∞),providedjk=1ukvkrevk=r!
,r=0,1,2j1.
(48)But(48)followsfrom(47),theequationsr!
=∞0xrexdx,r=0,1,∞0Lj(x)Lk(x)exdx=0forj=k,(see,forexample,[15,(6.
1.
1)and(22.
2.
13)]),andthefactthatGaussianquadratures,whichareexactforpolynomialsofdegreelessthantwicethenumberofnodes,havenodesthatcoincidewiththerootsofthecorrespondingorthogonalpolynomials(seesection2.
4).
Wehavecompletedthedenitionofthenewquadratures,alongwiththedemon-strationoftheirasymptoticperformance.
Weshallseethattheexistenceoftheserules,whichdependsonthesolvabilityofthenonlinearsystemsofequationsthatdenethenodesandweights,isassuredbythetheoryofChebyshevsystems.
Theuniquenessoftherulesissimilarlyassured.
Theseissuesofexistenceanduniquenessaretreatednext.
HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES15634.
Existenceanduniqueness.
4.
1.
Chebyshevsystems.
Materialofthissubsectionistaken,withminoralterations,fromKarlinandStudden[21].
SupposeIisanintervalofR,possiblyinnite.
Acollectionofnreal-valuedcontinuousfunctionsf1,fndenedonIisaChebyshevsystemifanylinearcombinationf(x)=ni=1aifi(x),withainotallzero,hasatmostn1zerosonI.
Thisconditionisequivalenttothestatementthatfordistinctx1,xninI,detf1(x1)···fn(x1).
.
.
.
.
.
.
.
.
f1(xn)···fn(xn)=0.
(49)TheChebyshevpropertyisacharacteristicofthespace,ratherthanthebasis:iff1,fnisaChebyshevsystem,thensoisanyotherbasisofspan{f1,fn}.
Ifuisacontinuous,positivefunctiononI,thenscalingbyupreservesaChebyshevsystem.
Finally,ifuisstrictlyincreasingandcontinuousonintervalJwithrangeI,thenf1u,fnuisaChebyshevsystemonJifandonlyiff1,fnisonI.
(Herefiudenotesthecompositionufollowedbyfi.
)Thebest-knownexampleofaChebyshevsystemisthesetofpolynomials1,x,xn1onanyintervalIR.
WeshallbeconcernedalsowiththeChebyshevsystems1,xγ,x,xγ+1x(n1)/2,xγ+(n1)/2,1,logx,x,xlogx,x(n1)/2,x(n1)/2logxonI=(0,a],whereγ∈R\Zanda>0.
ThesesystemsarespecialcasesofthesystemofM¨untzfunctions(see,forexample,[22,p.
133])M=xγilogkxk=0,ni1,i=1,j(50)onI=(0,∞),whereγ1,γjaredistinctrealnumbersandn1,njarepositiveintegerswithni=n.
ToseethisisaChebyshevsystem,supposef∈spanManduseinductioninnon(d/dlogx)[f(x)xγj],incombinationwithRolle'stheorem.
AnotherChebyshevsystemthatwillariseisthesystemL=(x+γi)kk=1,ni,i=1,j(51)onI=[0,∞),whereγ1,γjaredistinctpositiverealnumbersandn1,njarepositiveintegerswithni=n.
ThisisindeedaChebyshevsystem,foriff∈spanL,thenthefunctionf(x)ji=1(x+γi)niisapolynomialinxofdegreen1.
Supposef1,fnisaChebyshevsystemontheintervalI.
ThemomentspaceMnwithrespecttof1,fnisthesetMn=c=(c1,cn)∈Rnci=Ifi(x)dσ(x),i=1,n,(52)1564BRADLEYK.
ALPERTwherethemeasureσrangesoverthesetofnondecreasingright-continuousfunctionsofboundedvariationonI.
ItcanbeshownthatMnistheconvexconeassociatedwithpointsinthecurveC,whereC=f1(x)fn(x)x∈I.
Inotherwords,MncanberepresentedasMn=c=pj=1αjyjαj>0,yj∈C,j=1,p,p≥1.
TheindexI(c)ofapointcofMnistheminimumnumberofpointsofCthatcanbeusedintherepresentationofc,undertheconventionthatapoint(f1(x)fn(x))iscountedasahalfpointifxisfromtheboundaryofIandreceivesafullcountotherwise.
Theindexofaquadratureinvolvingx1,xpisdeterminedbycountinglikewise.
Proofsofthenextthreetheoremsaresomewhatelaborateandareomittedhere;theycanbefoundinKarlinandStudden[21].
Theorem4.
1.
(See[21,p.
42].
)SupposeI=[a,b]isaclosedinterval.
Apointc∈Mn,c=0,isaboundarypointofMnifandonlyifI(c)0,i=1,p.
Theorem4.
2.
(See[21,p.
47].
)SupposeI=[a,b]isaclosedinterval.
AnypointcintheinteriorofMnsatisesI(c)=n/2.
Moreover,ifσisameasurecorrespondingtoc,thenthereareexactlytwoquadraturespi=1wifr(xi)=Ifr(x)dσ(x),r=1,n,(54)ofindexn/2,wherewi>0,i=1,p.
Inparticular,ifn=2m,thenp=morp=m+1anda0implyni=1aici>0.
(61)Theorem4.
4.
(See[21,p.
106].
)Supposefi(x)=xi1fori=1,n,andI=[a,b].
Ifn=2m,thenc=(c1,cn)isanelementofMnifandonlyifthetwoquadraticformsmi,j=1[ci+jaci+j1]αiαjandmi,j=1[bci+j1ci+j]βiβj(62)arenonnegativedenite.
Ifn=2m+1,thenc∈Mnifandonlyifthetwoquadraticformsm+1i,j=1ci+j1αiαjandmi,j=1[(a+b)ci+jabci+j1ci+j+1]βiβj(63)arenonnegativedenite.
Moreover,foreitherparityofn,cisintheinteriorofMnifandonlyifthecorrespondingquadraticformsarebothpositivedenite.
Proof.
AtheoremofLukacs(see,forexample,[19,p.
4])statesthatapolynomialfofdegreen1thatisnonnegativeon[a,b]canberepresentedintheformf(x)=(xa)p(x)2+(bx)q(x)2,n=2m,p(x)2+(bx)(xa)q(x)2,n=2m+1,(64)wherepandqarepolynomialssuchthatthedegreeofeachtermin(64)doesnotexceedn1.
Thecombinationof(64)andTheorem4.
3provesthetheorem.
4.
2.
M¨untzsystemquadratures.
Thesystemsof(22),(24),(36),and(39)thatdenethequadraturerulesofsection3arespecialcasesofthesystemofequations(65)n/2m=1wmxmγilogkxm=(1)k+1ζ(k)(γi,a),k=0,ni1,i=1,j,fordistinctrealnumbersγ1,γjandpositiveintegersn1,njwithni=n.
Hereζ(k)denotesthekthderivativeofζwithrespecttoitsrstargument.
Theexis-tenceanduniquenessofthesolutionof(65)followfromtheexistenceanduniquenessofquadraturesforChebyshevsystems,onceitisestablishedthatthereisameasureσawith(66)a0xγilogkxdσa(x)=(1)k+1ζ(k)(γi,a),k=0,ni1,i=1,j,1566BRADLEYK.
ALPERTinotherwords,thatthemomentspaceMnoftheChebyshevsystemofM¨untzfunc-tionsM=xγilogkxk=0,ni1,i=1,j(67)on(0,a]containsthepointc=ζ(γ1,a)1)njζ(nj1)(γj,a).
(68)Wewillshowthatthisconditionissatisedprovidedthataissucientlylarge.
Itwouldbeconvenienttohavetightboundsfora,inparticularforsystems(22),(24),(36),and(39),butitappearsthatsuchboundsarediculttoobtain.
Evenfortheregularcases(22)and(24),wherebyTheorem4.
4theexistenceofσaisequivalenttothepositivedenitenessoftwomatrices,preciseboundsforarbitraryjappeardicult.
(Numericalexamplesbelowprovideevidencethata/jmaybechosenassmallas5/6.
)Theorem4.
5.
Supposeγ1,γjaredistinctrealnumbers,eachgreaterthan1,andn1,njarepositiveintegerswithni=n.
Thenforsucientlylargea,thereexistsameasureσasuchthatthesystemof(66)issatisedandcdenedby(68)isintheinteriorofthemomentspaceMn.
Proof.
Weconstructacontinuousweightfunctionσasatisfying(66)andshowthatforsucientlylargea,σa(x)ispositiveforx∈[0,a].
Welinearlycombinetheequationsof(66)toobtaintheequivalentsystem(69)a0(x/a)γilogk(x/a)dσa(x)=(1)k+1kr=0krζ(r)(γi,a)aγilogkra,k=0,ni1,i=1,j,wherewehaveusedthebinomialtheoremtoexpandlogk(x/a)=(logxloga)k.
Wedenetheweightσabytheformulaσa(x)=n1m=0αm,a(x/a)m(70)andcombine(69),(70),andtheequalities10xγlogkxdx=(1)kk!
(1+γ)k+1,γ>1,k=0,1,toobtaintheequationsinα0,a,αn1,a,(71)n1m=0αm,a(1)kk!
(1+γi+m)k+1=(1)k+1kr=0krζ(r)(γi,a)a1+γilogkra,k=0,ni1,i=1,j.
Thisn-dimensionallinearsystemisnonsingular,sincethesetoffunctions(x+γi+1)kk=1,ni,i=1,jformsaChebyshevsystemon[0,∞),asestablishedat(51).
Thus(71)possessesauniquesolutionα0,a,αn1,a.
HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1567Wenowdetermineαm=lima→∞αm,a,m=0,n1.
(72)Theasymptoticexpansionofζ(r)(γi,a)asa→∞canbederivedbydierentiating(10);therstseveraltermsaregivenbyζ(r)(γi,a)=a1+γirl=0rl(1)l+1(rl)!
logla(1+γi)rl+1+O(aγilogra).
(73)Combining(71)and(73),changingtheorderofsummation,andtwiceapplyingtheproductdierentiationruledrdγrf(γ)g(γ)=rs=0rsf(s)(γ)g(rs)(γ),weobtainnm=0αm(1)kk!
(1+γi+m)k+1=(1)kk!
(1+γi)k+1,k=0,ni1,i=1,j,whichimmediatelyreducestoαm=1,m=0,0,m=1,n1.
(74)Thecombinationof(70),(72),and(74)giveslima→∞σa(ax)=1,x∈[0,1],whichimpliesthatforasucientlylarge,σa(x)>0forx∈[0,a].
Thepointcdenedby(68)isintheinteriorofMn,sincesmallperturbationsofcwillpreservethepositivityofσa.
Theorem4.
2ensurestheexistenceofGaussianquadraturesforaChebyshevsys-temf1,fndenedonanintervalI,undertheassumptionthatIisclosed,whereasthesystemMof(67)isChebyshevonI=(0,a].
Asaconsequence,werequirethefollowingresult.
Theorem4.
6.
Supposethecollectionoffunctionsf1,fnformsaChebyshevsystemonI=(a,b]andeachisintegrableon[a,b]withrespecttoameasureσcorre-spondingtoapointcintheinteriorofMn.
Thenthereexistsauniquequadraturepi=1wifr(xi)=bafr(x)dσ(x),r=1,n,(75)ofindexn/2,wherewi>0andxi∈I,i=1,p.
Inparticular,ifn=2m,thenp=manda0fori=1,n/2and00andα>0suchthatthesequencex1,x2,.
.
.
generatedby(86)convergestoxandxi+1x≤αxix2(87)foranyinitialpointx0suchthatx0x0thereisaprocedurePtocomputextfort∈[0,1],givenanestimatextwith|xtxt|<δ.
ThenthereexistsapositiveintegermsuchthatthefollowingprocedurecanbeusedtocomputethesolutionofF(x)=0:Fori=1,m,usePtocomputexi/m,giventheestimatex(i1)/m.
TherequiredsolutionofF(x)=0isx1.
Moretypically,ofcourse,δandanyboundon|xt+xt|/dependontandinapracticalimplementationthestepsizeischosenadaptively.
Tocomputethesolutionsof(36)and(39),itiseectivetouseacontinuationprocedurewithrespecttobothjanda.
Solutionsfortherstfewvaluesofjarereadilyobtainedwithoutrequiringgoodinitialestimates.
Givenasolutionof(36)fortheinterval[0,a]withnodesu1,ujandweightsv1,vj,wechooseaninitialestimateu1,uj+1,v1,vj+1forj+1andtheinterval[0,a+1]denedbytheformulaeui=ui,i=1,j,a,i=j+1,vi=vi,i=1,j,1,i=j+1.
(89)Thischoiceexactlysatisestheequationsj+1i=1uiviγ+r=ζ(γr,a+1),r=0,1,j1,j+1i=1uivir=Br+1(a+1)r+1,r=0,1,j1,(90)asfollowsimmediatelyfromthedierenceformula(4)forBnand(11)forζ,butthecorrespondingequationsforr=jarenotsatised.
Thoseequationsareapproxi-matelysatised,however,andwecanstartwiththeactualvaluesofthesumsforr=jastherequiredvalues.
Thesearethenvariedcontinuously,obtainingthecor-respondingsolutions,untiltheycoincidewiththeintendedvaluesζ(γj,a+1)andBj+1(a+1)/(j+1).
Thisprocedurecanbeusedwithoutalterationfor(39).
Oncethesolutionforj+1andtheinterval[0,a+1]isobtained,acanbecontinuouslyvaried,inacontinuationprocedure,toobtainsolutionsfordierentintervals.
1572BRADLEYK.
ALPERTTable1Theminimumvalueofa,asafunctionofj,suchthatthepointB1(a)/1,Bj(a)/jisinthemomentspaceM2jofthepolynomials1,x,x2j1ontheinterval[0,a].
Themomentspaceisdenedat(52).
jminajminajminajmina10.
7886853.
9669697.
210811310.
4788521.
5708564.
77448108.
026181411.
2981532.
3634775.
58463118.
842741512.
1181543.
1628886.
39687129.
660351612.
938786.
Numericalexamples.
Theproceduresdescribedinsection5wereimple-mentedinPari/GP[23]forboththeregularcasesandthesingularcases.
Thematrixin(79),whichmustbeinverted,isverypoorlyconditionedformanychoicesofn,x1,xn,andf1,f2n.
ThisdicultywasmetbyusingtheextendedprecisioncapabilityofPari/GP.
6.
1.
Nodesandweights.
Thenodesandweightsof(22),(24),(36),and(39)thatdeterminethequadraturesofsection3werecomputedforarangeofvaluesoftheparameterj.
Foreachchoiceofj,awaschosen,byexperiment,tobethesmallestinte-gerleadingtopositivenodesandweights(seeTheorem4.
7).
Fortheregularcase(22),thecharacterizationexpressedinTheorem4.
4wasusedtodeterminetheminimumvalueofa∈R,forj=1,16,thatsatisesc=B1(a)/1,Bj(a)/j∈M2j.
Inparticular,weobtainedtheminimumvalueofasuchthatthequadraticformsin(62)arenonnegativedenite.
Thisdeterminationwasmadebycalculatingthede-terminantofeachcorrespondingmatrixsymbolicallyandsolvingforthelargestrootoftheresultingpolynomialina.
ThesevaluesaregiveninTable1.
Thisevidencesuggeststhatlimj→∞j1minaexistsandisroughly5/6,meaningthatthenumberoftrapezoidalnodesdisplacedislessthanthenumberofGaussiannodesreplacingthem,forquadraturerulesofallorders.
Thisrelationshipalsoappearstohold,toanevengreaterextent,forthesingularcases.
Thevaluesofselectednodesandweights,fortheregularcasesandforsingularitiesx1/2andlogx,arepresentedinanappendix.
Ofparticularsimplicityarethersttworulesforregularintegrands,T1,1n(f)=h12f(h/6)+f(h)f(1h)+12f(1h/6),(91)whereh=1/(n+1),forn=0,1,and(92)T2,2n(f)=h2548f(h/5)+4748f(h)+f(2h)+f(12h)+4748f(1h)+2548f(1h/5),whereh=1/(n+3),forn=0,1,Theserulesareofthird-andfourth-orderconvergence,respectively.
Therstisnoteworthyforhavingthesameweightsas,buthigherorderthan,thetrapezoidalrule;thesecondhasasymptoticerror1/4thatofSimpson'srulewiththesamenumberofnodes.
Thelowest-orderrulepresentedforlogarithmicsingularities,S1,1,1,1n(g)=h12g(h/(2π))+g(h)g(1h)+12g(1),(93)HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1573Table2Relativeerrorsinthecomputationoftheintegralin(94),fortheregularcases(x)=0.
Quadra-tureruleswithconvergenceoforder2,4,8,16,and32wereusedwithvariousnumbersm=n+2jofnodes.
Heref=m/(200/π)istheoversamplingfactor.
mf2481632701.
100.
622D+000.
114D010.
382D020.
170D050.
234D10801.
260.
488D+000.
938D020.
184D020.
354D060.
115D11901.
410.
391D+000.
744D020.
934D030.
841D070.
720D131001.
570.
321D+000.
584D020.
498D030.
223D070.
192D141151.
810.
246D+000.
408D020.
211D030.
365D080.
331D141302.
040.
194D+000.
289D020.
964D040.
715D090.
331D141452.
280.
157D+000.
209D020.
472D040.
162D090.
471D141602.
510.
129D+000.
154D020.
245D040.
415D100.
262D141802.
830.
102D+000.
106D020.
110D040.
794D110.
471D142003.
140.
832D010.
747D030.
531D050.
177D110.
331D142303.
610.
631D010.
465D030.
199D050.
235D120.
523D152604.
080.
495D010.
303D030.
828D060.
375D130.
384D14Table3Relativeerrorsforthesingularcases(x)=x1/2,forvariousnumbersm=n+j+kofnodesandordersofconvergence.
mf24816701.
100.
692D010.
519D010.
850D020.
163D03801.
260.
925D010.
258D010.
260D020.
578D05901.
410.
921D010.
133D010.
698D030.
667D061001.
570.
838D010.
717D020.
146D030.
277D061151.
810.
686D010.
307D020.
201D040.
360D071302.
040.
550D010.
144D020.
269D040.
437D081452.
280.
441D010.
724D030.
171D040.
557D091602.
510.
357D010.
389D030.
964D050.
733D101802.
830.
273D010.
186D030.
440D050.
408D112003.
140.
212D010.
976D040.
207D050.
218D122303.
610.
151D010.
427D040.
724D060.
130D122604.
080.
110D010.
215D040.
280D060.
201D13approximates10g(x)dxwitherroroforderO(h2logh)forg(x)=φ(x)logx+ψ(x),providedφandψareregularfunctionson[0,1].
Thecorrespondingruleforthesingularityx1/2isnotquiteassimple,fortherethequantity2πin(93)isreplacedwith4ζ(1/2)2.
6.
2.
Quadratureperformance.
Todemonstratetheperformanceofthequad-raturerules,theywereusedinaFortranroutine(withreal*8arithmetic)tonumer-icallycomputetheintegrals10cos(200x)s(x)+cos(200x+.
3)dx,(94)forthefunctionss(x)=0,s(x)=x1/2,ands(x)=logx.
Theseintegralswerealsoobtainedanalyticallyandtherelativeerrorofthequadratureswascomputed.
Thenumericalintegrationswerecomputedforvariousordersofquadratureandvariousnumbersofnodes.
Minimumsamplingwastakentobetwopointsperperiodofthecosine(i.
e.
,200/π≈63.
7quadraturenodes).
Theaccuracieswerethencomparedforvariousdegreesofoversampling.
ThequadratureerrorsarelistedinTables2–4andplotted,asafunctionofoversamplingfactor,inFigure1.
Wenotethatthegraphsare1574BRADLEYK.
ALPERTTable4Relativeerrorsforthesingularcases(x)=logx.
HeretheerrorisoforderO(hllogh),wherelisshown.
mf24816701.
100.
369D+000.
217D010.
354D010.
243D03801.
260.
271D+000.
238D020.
328D020.
487D04901.
410.
206D+000.
765D020.
707D030.
394D051001.
570.
162D+000.
768D020.
687D030.
121D051151.
810.
117D+000.
576D020.
291D030.
886D071302.
040.
882D010.
398D020.
120D030.
903D081452.
280.
687D010.
272D020.
548D040.
123D081602.
510.
549D010.
188D020.
272D040.
177D091802.
830.
421D010.
119D020.
118D040.
965D112003.
140.
332D010.
774D030.
550D050.
956D122303.
610.
243D010.
433D030.
196D050.
398D122604.
080.
185D010.
258D030.
778D060.
106D1212341016101410121010108106104102100s(x)=0Oversampling(f)RelativeError24816321234s(x)=x^(1/2)Oversampling(f)248161234s(x)=log(x)Oversampling(f)24816Fig.
1.
Therelativeerrorsofthequadratures,showninTables2–4,areplottedusinglogarithmicscalingonbothaxes.
nearlystraightlines(untilthelimitofmachineprecisionisreached),aspredictedfromthetheoreticalconvergencerates.
Weremarkalsothatexcellentaccuracyisattainedforevenquitemodestoversamplingwhenquadratureswithhigh-orderconvergenceareemployed.
Forproblemswherethenumberofquadraturenodesisthemajorcostfactor,therefore,onemaybenetbyusingthehigh-orderquadraturesevenformodestaccuracyrequirements.
HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1575Table5Relativeerrorsinthecomputationoftheintegralin(95),forξ=1.
QuadraturerulesdenedinCorollary3.
10withj=1,2,4,8,and16wereusedwithvariousnumbersmofnodes.
m124816700.
999D+000.
400D+000.
305D+000.
180D+000.
104D01800.
304D+000.
200D010.
247D010.
238D020.
474D03900.
113D+000.
217D010.
136D020.
383D030.
866D051000.
273D010.
137D010.
137D020.
440D040.
900D061150.
210D010.
247D020.
137D030.
573D050.
331D071300.
228D010.
107D020.
632D040.
423D060.
107D091450.
118D010.
115D020.
305D040.
196D060.
490D101600.
212D020.
521D030.
307D050.
430D070.
216D091800.
625D020.
824D040.
451D050.
101D070.
166D112000.
558D020.
206D030.
184D050.
291D080.
867D122300.
863D030.
676D040.
463D060.
626D090.
586D132600.
266D020.
433D040.
291D060.
596D100.
346D14Wetestthequadraturesforimproperintegralsbynumericallycomputingforξ=1theintegral∞∞eixξ10r=10r+1x+r+idx=2πiH(ξ)10r=10(r+1)eirξξ,(95)whereHistheHeavisidestepfunction.
Theintegrandisoscillatoryanddecayslikex1asx→±∞.
ThequadraturesdenedinCorollary3.
10areemployed,forwhichtheintegralissplitintoaregularintegralonaniteinterval,chosenheretobe[5√m/4,5√m/4],wheremisthetotalnumberofquadraturenodes,andtwoimproperintegralsintheimaginarydirection,usingLaguerrequadratures.
ThequadratureerrorsareshowninTable5.
7.
Applicationsandsummary.
ThechiefmotivationforthehybridGauss-trapezoidalquadraturerulesistheaccuratecomputationofintegraloperators.
WedeneanintegraloperatorAbytheformula(Af)(x)=ΓK(x,y)f(y)dy,whereΓisaregular,simpleclosedcurveinthecomplexplane,thefunctionfisregularonΓ,andthekernelK:C*C→Cisaregularfunctionofitsarguments,exceptwheretheycoincide;weassumeK(x,y)=φ(x,y)s(|xy|)+ψ(x,y)(96)withφandψregularonΓ*Γandsregularon(0,∞),withanintegrablesingularityat0.
Alargevarietyofproblemsofclassicalphysicscanbeformulatedasintegralequationsthatinvolvesuchoperators.
Whentheoperatoroccursinanintegralequa-tionf(x)+(Af)(x)=g(x),x∈Γ,(97)somechoiceofdiscretizationmustbeusedtoreducetheproblemtoanite-dimensionalonefornumericalsolution.
IntheNystr¨ommethodtheintegralsarereplacedbyquadraturestoyieldthenitesystemofequationsf(xi)+mj=1wijf(xj)=g(xi),i=1,m.
(98)1576BRADLEYK.
ALPERTThislinearsystemcanbesolvedforf(x1)f(xm)byavarietyoftechniques.
Theparticularchoiceofxiandwijfori,j=1,mdeterminestheorderofconvergence(andthereforeeciency)ofthemethod.
Foracurveparametrizationν:[0,1]→Γ,suchasscaledarclength,theoperatorAbecomes(Af)(ν(t))=10K(ν(t),ν(τ))f(ν(τ))ν(τ)dτ.
Itisconvenienttouseauniformdiscretization1/m,2/m,1intandτ,soxi=ν(i/m),i=1,m.
HowtheniswijdeterminedWeassumeforthemomentthatfisavailableatlocationsotherthanx1,xm.
Continuingνperiodicallywithperiod1,andusingtheGauss-trapezoidalquadratures,weobtain(99)(Af)(ν(i/m))=1+i/mi/mK(ν(i/m),ν(τ))f(ν(τ)ν(τ))dτ≈1mjk=1ukσi/m(vk/m)+1mn1k=0σi/m(a/m+k/m)+1mjk=1ukσi/m(1vk/m),fori=1,m,whereσα:[0,1]→Cisdenedbytheformulaσα(τ)=K(ν(α),ν(α+τ))f(α+τ)ν(α+τ)(100)andm=n+2a1andu1,uj,v1,vjaredeterminedforthesingularitysofK.
Providedthattheperiodiccontinuationofνissucientlyregular,thequadraturewillconvergetotheintegralwithordergreaterthanjasm→∞,fori=1,m.
Werelaxtherestrictionthatfbeavailableoutsidex1,xmbyusinglocalLagrangeinterpolationoforderj+1forequispacednodes,f(ν(τ))≈jr=0f(ν(i/m+r/m))lr(mτi),(101)wherei=mτ(j1)/2andlr(x)=js=0,s=rxsrs,r=0,j.
Nowwijfori,j=1,misdeterminedbycombining(97)–(101).
Thecomputationofallm2coecientsrequiresm(m+2j2a+1)evaluationsofthekernelKandthereforeorderO(m2)operations.
Thiscostcanoftenbesubstantiallyreducedusingtechniquesthatexploitkernelsmoothness(see,forexample,[6],[5]).
AslightlydierentapplicationofthequadraturesisthecomputationofFouriertransformsoffunctionsthatfailtosatisfytheassumptionsusuallymadewhenusingthediscreteFouriertransform.
Inparticular,ifafunctiondecaysslowlyforlargeargumentoriscompactlysupportedandsingularattheendsofthesupportinterval,thesequadraturescanbeusedtocomputeitsFouriertransform.
Oneexampleofsuchafunctionisthatin(95).
Sincemostofthenodesinthesequadraturesareequispaced,withfunctionvaluesgivenequalweight,thefastFouriertransformcanbeusedtodothebulkofthecomputations;theoverallcomplexityisO(nlogn),wherenisthenumberofFouriercoecientstobecomputed.
OtherapplicationsmayincludetheHYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1577representationoffunctionsforsolvingordinaryorpartialdierentialequations,whenhigh-ordermethodsarerequired.
Inaddition,anextensionofthesequadraturestointegralsonsurfacesisunderstudy.
Insummary,thecharacteristicsofthehybridGauss-trapezoidalquadraturerulesincludearbitraryorderconvergenceforregularfunctionsorfunctionswithknownsingularitiesofpowerorlogarithmictype,positiveweights,mostnodesequispacedandmostweightsconstant,andinvariantnodesandweights(asidefromscaling)astheproblemsizeincreases.
Theprimarydisadvantageofthequadraturerules,sharedwithotherGaussianquadra-turesbutexacerbatedherebypoorconditioning,isthatthecomputationofthenodesandweightsisnottrivial.
Nevertheless,tabulationofnodesandweightsforagivenorderofconvergenceallowsthisissuetobeavoidedintheconstructionofhigh-order,general-purposequadratureroutines.
Appendix.
Tablesofquadraturenodesandweights.
Quadraturenodesandweightsmayalsobeobtainedelectronicallyfromtheau-thor.
1578BRADLEYK.
ALPERTTable6ThenodesandweightsforthequadratureruleTjan(f)=hji=1wif(xih)+hn1i=0f(ah+ih)+hji=1wif(1xih),withh=(n+2a1)1,forseveralchoicesofjandcorrespondingminimumintegera.
Forfaregularfunction,Tjan(f)convergesto10f(x)dxasn→∞withconvergenceoforderO.
Oaxiwi311.
666666666666667D015.
000000000000000D01422.
000000000000000D015.
208333333333333D011.
000000000000000D+009.
791666666666667D01522.
245784979812614D015.
540781643606372D011.
013719374359164D+009.
459218356393628D01632.
250991042610971D015.
549724327164180D011.
014269060987992D+009.
451317411845473D012.
000000000000000D+009.
998958260990347D01732.
180540672543505D015.
408088967208193D011.
001181873031216D+009.
516615045823566D011.
997580526418033D+001.
007529598696824D+00842.
087647422032129D015.
207988277246498D019.
786087373714483D019.
535038018555888D011.
989541386579751D+001.
024871626402471D+003.
000000000000000D+001.
000825744017291D+001257.
023955461621939D021.
922315977843698D014.
312297857227970D015.
348399530514687D011.
117752734518115D+008.
170209442488760D012.
017343724572518D+009.
592111521445966D013.
000837842847590D+009.
967143408044999D014.
000000000000000D+009.
999820119661890D011679.
919337841451028D022.
528198928766921D015.
076592669645529D015.
550158230159486D011.
184972925827278D+007.
852321453615224D012.
047493467134072D+009.
245915673876714D013.
007168911869310D+009.
839350200445296D014.
000474996776184D+009.
984463448413151D015.
000007879022339D+009.
999592378464547D016.
000000000000000D+009.
999999686258662D012099.
209200446233291D022.
351836144643984D014.
752021947758861D015.
248820509085946D011.
124687945844539D+007.
634026409869887D011.
977387385642367D+009.
284711336658351D012.
953848957822108D+001.
010969886587741D+003.
976136786048776D+001.
024959725311073D+004.
994354281979877D+001.
010517534639652D+005.
999469539335291D+001.
001551595797932D+006.
999986704874333D+001.
000061681794188D+008.
000000000000000D+001.
000000135843597D+00HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1579Table6(Continued)Oaxiwi24106.
001064731474805D021.
538932104518340D013.
149685016229433D013.
551058128559424D017.
664508240518316D015.
449200036280007D011.
396685781342510D+007.
104078497715549D012.
175195903206602D+008.
398780940253654D013.
062320575880355D+009.
272767950890611D014.
016440988792476D+009.
750605697371132D015.
002872064275734D+009.
942629650823470D016.
000285453310164D+009.
992421778421898D017.
000012964962529D+009.
999534370786161D018.
000000175554469D+009.
999990854912925D019.
000000000000000D+009.
999999989466828D0128126.
234360533194102D021.
595975279734157D013.
250286721702614D013.
637046028193864D017.
837350794282182D015.
498753177297441D011.
415673112616924D+007.
087986792086956D012.
189894250061313D+008.
335172275501195D013.
070053877483040D+009.
204446510608518D014.
018613756218047D+009.
710881776552090D015.
002705902035397D+009.
933296578555239D015.
999929741810400D+009.
994759087910050D016.
999904720846024D+001.
000133030254421D+007.
999986894843540D+001.
000032915011460D+008.
999999373380393D+001.
000002261653775D+009.
999999992002911D+001.
000000042393520D+001.
100000000000000D+011.
000000000042872D+0032145.
899550614325259D021.
511076023874179D013.
082757062227814D013.
459395921169090D017.
463707253079130D015.
273502805146873D011.
355993726494664D+006.
878444094543021D012.
112943217346336D+008.
210319140034114D012.
987241496545946D+009.
218382875515803D013.
944798920961176D+009.
873027487553060D014.
950269202842798D+001.
018251913441155D+005.
972123043117706D+001.
021933430349293D+006.
989783558137742D+001.
012567983413513D+007.
997673019512965D+001.
004052289554521D+008.
999694932747039D+001.
000713413344501D+009.
999979225211805D+001.
000063618302950D+001.
099999938266130D+011.
000002486385216D+001.
199999999462073D+011.
000000030404477D+001.
300000000000000D+011.
000000000020760D+001580BRADLEYK.
ALPERTTable7Thenodesv1,vjandweightsu1,ujforthequadratureruleSjkabn(g)=hji=1uig(vih)+hn1i=0g(ah+ih)+hki=1wig(1xih),withh=(n+a+b1)1,forg(x)=x1/2φ(x)+ψ(x),withφandψregularfunctions.
Thenodesx1,xkandweightsw1,wkarefoundinTable6.
Oaviui1.
511.
172258571393266D015.
000000000000000D012.
029.
252112715421378D024.
198079625266162D011.
000000000000000D001.
080192037473384D+002.
526.
023873796408450D022.
858439990420468D018.
780704050676215D011.
214156000957953D+003.
027.
262978413470474D033.
907638767531813D022.
246325512521893D014.
873484056646474D011.
000000000000000D+009.
735752066600344D013.
521.
282368909458828D026.
363996663105925D022.
694286346792474D015.
077434578043636D011.
018414523786358D+009.
286165755645772D014.
031.
189242434021285D025.
927215035616424D022.
578220434738662D014.
955981740306228D011.
007750064585281D+009.
427131290628058D012.
000000000000000D+001.
002416546550407D+006.
043.
317925942699451D031.
681780929883469D028.
283019705296352D021.
755244404544475D014.
136094925726231D015.
039350503858001D011.
088744373688402D+008.
266241339680867D012.
006482101852379D+009.
773065848981277D013.
000000000000000D+009.
997919809947032D018.
051.
214130606523435D036.
199844884297793D033.
223952700027058D027.
106286791720044D021.
790935383649920D012.
408930104410471D015.
437663805244631D014.
975929263668960D011.
176116628396759D+007.
592446540441226D012.
031848210716014D+009.
322446399614420D013.
001961225690812D+009.
928171438160095D014.
000000000000000D+009.
999449125689846D0110.
061.
745862989163252D041.
016950985948944D038.
613670540457314D032.
294670686517670D026.
733385088703690D021.
076657968022888D012.
514488774733840D012.
734577662465576D016.
341845573737690D014.
978815591924992D011.
248404055083152D+007.
256208919565360D012.
065688031953401D+008.
952638690320078D013.
009199358662542D+009.
778157465381624D014.
000416269690208D+009.
983390781399277D015.
000000000000000D+009.
999916342408948D01HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1581Table7(Continued)Oaviui12.
085.
710218427206990D042.
921018926912141D031.
540424351115548D023.
431130611256885D028.
834248407196555D021.
224669495638615D012.
824462054509770D012.
761108242022520D016.
574869892305580D014.
797809643010337D011.
246541060977993D+006.
966555677271379D012.
039218495130811D+008.
790077941972658D012.
979333487049800D+009.
868622449294327D013.
985772595393049D+001.
015142389688201D+004.
997240804311428D+001.
006209712632210D+005.
999868793951190D+001.
000528829922287D+007.
000000000000000D+001.
000002397796838D+0014.
093.
419821460249725D041.
750957243202047D039.
296593430187960D032.
080726584287380D025.
406214771755252D027.
586830616433430D021.
763945096508648D011.
766020526671851D014.
218486605653738D013.
206624362072232D018.
274022895884040D014.
934405290553812D011.
410287585637014D+006.
707497030698472D012.
160997505238153D+008.
244959025366557D013.
043504749358223D+009.
314646742162802D014.
005692579069439D+009.
845768443163154D014.
999732707905968D+009.
992852769154770D015.
999875191971098D+001.
000273112957723D+006.
999994560568667D+001.
000022857402321D+008.
000000000000000D+001.
000000081405180D+0016.
0102.
158438988280793D041.
105804873501181D035.
898432743709196D031.
324499944707956D023.
462795956896131D024.
899842307592144D021.
145586495070213D011.
165326192868815D012.
790344218856415D012.
178586693194957D015.
600113798653321D013.
481766016945031D019.
814091242883119D014.
964027915911545D011.
553594853974655D+006.
469026189623831D012.
270179114036658D+007.
823688971783889D013.
108234601715371D+008.
877772445893361D014.
032930893996553D+009.
551665077035583D015.
006803270228157D+009.
876285579741800D016.
000815466735179D+009.
979929183863017D017.
000045035079542D+009.
998470620634641D018.
000000738923901D+009.
999962891645340D019.
000000000000000D+009.
999999946893169D011582BRADLEYK.
ALPERTTable8Thenodesv1,vjandweightsu1,ujforthequadratureruleSjkabn(g)=hji=1uig(vih)+hn1i=0g(ah+ih)+hki=1wig(1xih),withh=(n+a+b1)1,forg(x)=φ(x)logx+ψ(x),withφandψregularfunctions.
TheerrorisoforderO(hllogh).
Thenodesx1,xkandweightsw1,wkarefoundinTable6.
laviui211.
591549430918953D015.
000000000000000D01321.
150395811972836D013.
913373788753340D019.
365464527949632D011.
108662621124666D+00422.
379647284118974D028.
795942675593887D022.
935370741501914D014.
989017152913699D011.
023715124251890D+009.
131388579526912D01532.
339013027203800D028.
609736556158105D022.
854764931311984D014.
847019685417959D011.
005403327220700D+009.
152988869123725D011.
994970303994294D+001.
013901778984250D+00634.
004884194926570D031.
671879691147102D027.
745655373336686D021.
636958371447360D013.
972849993523248D014.
981856569770637D011.
075673352915104D+008.
372266245578912D012.
003796927111872D+009.
841730844088381D01856.
531815708567918D032.
462194198995203D029.
086744584657729D021.
701315866854178D013.
967966533375878D014.
609256358650077D011.
027856640525646D+007.
947291148621894D011.
945288592909266D+001.
008710414337933D+002.
980147933889640D+001.
036093649726216D+003.
998861349951123D+001.
004787656533285D+001061.
175089381227308D034.
560746882084207D031.
877034129831289D023.
810606322384757D029.
686468391426860D021.
293864997289512D013.
004818668002884D012.
884360381408835D016.
901331557173356D014.
958111914344961D011.
293695738083659D+007.
077154600594529D012.
090187729798780D+008.
741924365285083D013.
016719313149212D+009.
661361986515218D014.
001369747872486D+009.
957887866078700D015.
000025661793423D+009.
998665787423845D01HYBRIDGAUSS-TRAPEZOIDALQUADRATURERULES1583Table8(Continued)laviui1271.
674223682668368D036.
364190780720557D032.
441110095009738D024.
723964143287529D021.
153851297429517D011.
450891158385963D013.
345898490480388D013.
021659470785897D017.
329740531807683D014.
984270739715340D011.
332305048525433D+006.
971213795176096D012.
114358752325948D+008.
577295622757315D013.
026084549655318D+009.
544136554351155D014.
003166301292590D+009.
919938052776484D015.
000141170055870D+009.
994621875822987D016.
000001002441859D+009.
999934408092805D011499.
305182368545380D043.
545060644780164D031.
373832458434617D022.
681514031576498D026.
630752760779359D028.
504092035093420D021.
979971397622003D011.
854526216643691D014.
504313503816532D013.
251724374883192D018.
571888631101634D014.
911553747260108D011.
434505229617112D+006.
622933417369036D012.
175177834137754D+008.
137254578840510D013.
047955068386372D+009.
235595514944174D014.
004974906813428D+009.
821609923744658D014.
998525901820967D+001.
000047394596121D+005.
999523015116678D+001.
000909336693954D+006.
999963617883990D+001.
000119534283784D+007.
999999488130134D+001.
000002835746089D+0016108.
371529832014113D043.
190919086626234D031.
239382725542637D022.
423621380426338D026.
009290785739468D027.
740135521653088D021.
805991249601928D011.
704889420286369D014.
142832599028031D013.
029123478511309D017.
964747731112430D014.
652220834914617D011.
348993882467059D+006.
401489637096768D012.
073471660264395D+008.
051212946181061D012.
947904939031494D+009.
362411945698647D013.
928129252248612D+001.
014359775369075D+004.
957203086563112D+001.
035167721053657D+005.
986360113977494D+001.
020308624984610D+006.
997957704791519D+001.
004798397441514D+007.
999888757524622D+001.
000395017352309D+008.
999998754306120D+001.
000007149422537D+001584BRADLEYK.
ALPERTAcknowledgments.
TheauthorthanksGregoryBeylkin,VladimirRokhlin,andRonaldWittmannforhelpfuldiscussionsandencouragement.
REFERENCES[1]V.
Rokhlin,Rapidsolutionofintegralequationsofclassicalpotentialtheory,J.
Comput.
Phys.
,60(1983),pp.
187–207.
[2]V.
Rokhlin,DiagonalformsoftranslationoperatorsfortheHelmholtzequationinthreedi-mensions,Appl.
Comput.
Harmon.
Anal.
,1(1993),pp.
82–93.
[3]L.
GreengardandV.
Rokhlin,Afastalgorithmforparticlesimulations,J.
Comput.
Phys.
,73(1987),pp.
325–348.
[4]W.
HackbuschandZ.
Nowak,Onthefastmatrixmultiplicationintheboundaryelementmethodbypanelclustering,Numer.
Math.
,54(1989),pp.
463–491.
[5]G.
Beylkin,R.
Coifman,andV.
Rokhlin,Fastwavelettransformsandnumericalalgo-rithmsI,Comm.
PureAppl.
Math.
,XLIV(1991),pp.
141–183.
[6]B.
Alpert,G.
Beylkin,R.
Coifman,andV.
Rokhlin,Wavelet-likebasesforthefastsolutionofsecond-kindintegralequations,SIAMJ.
Sci.
Comput.
,14(1993),pp.
159–184.
[7]C.
T.
Kelley,Afastmultilevelalgorithmforintegralequations,SIAMJ.
Numer.
Anal.
,32(1995),pp.
501–513.
[8]I.
Navot,AnextensionoftheEuler-Maclaurinsummationformulatofunctionswithabranchsingularity,J.
Math.
Phys.
,40(1961),pp.
271–276.
[9]V.
Rokhlin,End-pointcorrectedtrapezoidalquadraturerulesforsingularfunctions,Comput.
Math.
Appl.
,20(1990),pp.
51–62.
[10]H.
P.
Starr,OntheNumericalSolutionofOne-DimensionalIntegralandDierentialEqua-tions,Ph.
D.
thesis,YaleUniversity,NewHaven,CT,1991.
[11]B.
Alpert,High-orderquadraturesforintegraloperatorswithsingularkernels,J.
Comput.
Appl.
Math.
,60(1995),pp.
367–378.
[12]R.
Kress,ANystr¨ommethodforboundaryintegralequationsindomainswithcorners,Numer.
Math.
,58(1990),pp.
145–161.
[13]S.
KapurandV.
Rokhlin,High-ordercorrectedtrapezoidalrulesforsingularfunctions,SIAMJ.
Numer.
Anal.
,34(1997),pp.
1331–1356.
[14]J.
Ma,V.
Rokhlin,andS.
Wandzura,GeneralizedGaussianquadraturerulesforsystemsofarbitraryfunctions,SIAMJ.
Numer.
Anal.
,33(1996),pp.
971–996.
[15]M.
AbramowitzandI.
A.
Stegun,Eds.
,HandbookofMathematicalFunctions(10thprint-ing),AppliedMathematicsSeries55,NationalBureauofStandards,1972.
[16]A.
Erdelyi,Ed.
,HigherTranscendentalFunctions,Vol.
1,Krieger,Malabar,FL,1985.
[17]C.
M.
BenderandS.
A.
Orszag,AdvancedMathematicalMethodsforScientistsandEngi-neers,McGraw-Hill,NewYork,1978.
[18]J.
StoerandR.
Bulirsch,IntroductiontoNumericalAnalysis,Springer-Verlag,Berlin,1980.
[19]G.
Szeg¨o,OrthogonalPolynomials,AmericanMathematicalSociety,Providence,RI,1939.
[20]J.
N.
Lyness,Finite-partintegralsandtheEuler-Maclaurinexpansion,inApproximationandComputation,Internat.
Ser.
Numer.
Math.
119,R.
V.
M.
Zahar,ed.
,Birkh¨auserVerlag,Basel,Boston,Berlin,1994,pp.
397–407.
[21]S.
KarlinandW.
Studden,TchebycheSystemswithApplicationsinAnalysisandStatistics,JohnWiley,NewYork,1966.
[22]P.
BorweinandT.
Erdelyi,PolynomialsandPolynomialInequalities,Springer-Verlag,NewYork,1995.
[23]C.
Batut,D.
Bernardi,H.
Cohen,andM.
Olivier,User'sGuidetoPARI-GP,Technicalre-port,LaboratoireA2X,UniversiteBordeaux,I,1995.
Availableviaftpfrommegrez.
math.
u-bordeaux.
fr.
[24]W.
Gautschi,OntheconstructionofGaussianquadraturerulesfrommodiedmoments,Math.
Comput.
,24(1970),pp.
245–260.
[25]G.
H.
GolubandJ.
H.
Welsch,CalculationofGaussianquadraturerules,Math.
Comput.
,23(1969),pp.
221–230.