C. R. Acad. Sci. Paris, t. 327, Analyse num&iquelNumerical
Inthgration r&guli~res Philippe
Serie
I, p. 843-848,
numhrique ou singuli&res
HELLUY
1998
Analysis
a, Sylvain
MAIRE
d’ordre Aeva de fonctions sur un intervalle a, Patrick
RAVEL
’ Laboratoire modblisation numbrique et couplage, institut Var, B.P. 56, 83162 La Valette du Var cedcx, Franc? Courriel : helluy@isitv. univ-tln.fr,
[email protected] Ir Lahoratoire pharmacie, Courriel (Rqu
de physique industriellr 34060 Montpellier cedrx :
[email protected]
le 3 juillet
RCsumiS.
1998,
ac:erptC:
aprRs
et traitemcnt 1. France
r&vision
des
sciences
de l’information,
Ir 28 septamhre
’ de l’ing&ieur
CR&A
URA
CNRS
de Toulon
1465,
facultk
et du
de
1998)
Nous prksentons une mkthode d’intkgration numkrique d’ordre Clev6 pour des fonctions singulibres ou non sur un intervalle. Cette mkthode s’avk-e trbs utile dans la mise en ceuvre des mkthodes d’C1Cments finis front&e, pour intkgrer la fonction de Green. 0 AcadCmie des SciencesElsevier, Paris
High
order
functions
numerical integration on an interval
of regular
or singular
Abstract.
A high order numerical integration method is presented for regular or singular integrands over an interval. This method appears to be very useful to compute the integrals of the Green ,function in the numerical resolution of boundarq integral equations. 0 AcadCmie des Sciences/Elsevier, Paris
A bridged
English
Version
It is important, for the numerical resolution of boundary integral equations, to construct numerical integration methods for regular or singular functions on an interval. Generally, classical Gauss formula are used in the domain of regularity of the integrand, and a special treatment is performed near the singularities to preserve high order. This treatment can lead to very tedious computations: elimination of the singularity, changesof variables, computation of special Gauss points and weights which depend on the nature of the singularity... We present here a very simple numerical quadrature method which is of high order both for regular and singular integrands. This method is based on the fact that the trapezoidal rule (or the rectangular Note prCsent6e par Philippe 0764-4442/98/03270843
G. CIARLET.
0 Acadkmie
des Sciences/Elsevier.
Paris
843
P. Helluy
et al.
rule) is of high order for regular periodic functions. The high order can be preserved with the help of a change of variables even when the integrand is no more periodic. This fact has been first discovered by hi, Moriguti and Takasawa in [3]. It is also described in [ 11. For more references, see also [4], [5]. In this note we study the case of a polynomial change of variable. We perform error estimates both for regular or singular integrands, and we prove numerically that our estimates are optimal. Our presentation is organised as follow: We first recall some well-known results about the trapezoidal rule (see formula (I)) applied to periodic functions. Those results are summarized in Theorem 1 where the error estimate (3) is proved for functions verifying (2). This estimate is based on the important expansion (4). We then extend the method to non-periodic functions. First, a polynomial change of variables is constructed. The polynom has to verify (5) and (6) and is given by formula (7). Then the trapezoidal rule is applied to P(f(t))P’(t) (see (8)), and it defines what we call the periodisation method. We are able to show the error estimate (9). Finally, we extend our method to singular integrands of the form (10) and are able to prove estimate (12). When the order of the singularity increases the order of the periodisation method decreases and a higher order of the polynom P has to be used. But the order is still higher than one (case of the classical Gauss method applied to singular integrands). Our numerical results confirm the previous theory. This method will be used to solve boundary integral equations in a forthcoming paper of Helluy, Maire and Ravel [2].
1. Ordre
de la mkthode
des rectangles
Soit f une fonction definie sur [0, 11, de classe C1 et soit J.~(.lf> (N entier > 0) l’approximation par la methode des rectangles a gauche de I(f) = J% f(C)dt
(1)
En general, cette methode d’integration numerique est d’ordre 1, c’est-a-dire que l’erreur EN(~) = r(f) - Jdf) d Ccroit comme & quand N tend vers +oe. Cependant, si f est plus reguliere et si ses derivees d’ordre impair vbifient certaines proprietes, la methode est d’ordre plus Cleve. Nous rappelons ce resultat dans le theoreme suivant. TH~OR~ZME 1. - Soit k un entier 2 2! et soit f me ,fonction de clusse C2” sur l’intervalle On suppose que
f’2i-lyq
=
f(2i4)
(1) pour 1 < i 5 k - 1,
alors il existe une constunte CI, > 0. indt?pendunte de f, telle que
844
[0, 11.
(2)
lnt+yation
Dkmonstration. - Par developpement on peut montrer clue
numkique
de fonctions
singulikes
ou
non
en serie de Fourier de f et integrations par parties successives,
(4)
la condition
(2) implique
done que
d’oti le resultat.
0
Remarque. - La demonstration ci-dessus utilise la convergence f. Ceci explique que l’estimation soit fausse pour k = 1. 2. IntCgration
des fonctions
rCguli&res
normale de la serie de Fourier de
sur [O, l] par la mkthode
de pkriodisation
Dans ce paragraphe nous introduisons la methode dite de periodisation. Un changement de variable bien choisi va now permettre de nous ramener au cas precedent. m&me lorsque la fonction f de classe C*” (k > 2) sur [0, l] ne verifie plus les conditions (2) sur ses derivees d’ordre impair. Pour cela now commenqons par construire un polynome P verifiant les proprietes suivantes : P(0) P(“)(O)
= 0,
P(1) = 1.
= P(‘)(l)
= 0,
(5) 1 < 15 2k - 2.
(6)
LEMME 1. - I1 existe un unique polyncime P de degre’ 4k - 3 ve’rijiant (5) et (6). Ce polyncime est strictement croissant sur [0, 11. I1 est donnP par la formule P(X) = La methode de periodisation oti 9 est la fonction
.I;” P-2(
1 - t)2”-2dt
\’ t2’;-2(] .0
_ Q2”-2&’
consiste alors a remarquer
que 1(f)
(7) = Jo1f(P(u))P’(u)ciw,
= I(g),
.9(u) = fuw)~‘(4. Par construction de P, il se trouve que g et un certain nombre de ses derivees s’annulent en 0 1. 11 suffit done d’appliquer la methode des rectangles composite a 9 pour construire une methode d’integration d’ordre Clew? : et
LEMME
lafonction
2. - Soit f : [0, l] -+ [w unefonction de classe C *’ . Soit P le polyncime donne’ par (7). Alors 9 dkjinie sur [O: l] par 9(u) = f(P(u))P’(u) ve’rijie g(J)(O) = p(1)
= 0 pour 0 <-. j <- 2x: - 3
845
P. Helluy
et al.
COROLLAIRE 1. - Soit f : [O, l] -+ W we fonction de classe C 2k. Soit P le polynBme don& par (7).
On pose
alors il existe une constante GI, ne d&pendant que de k telle que
3. GCnhlisation
A des fonctions
singuliirres
aux bornes de I’intervalle
d’intbgration
Un des avantages de la methode de periodisation est qu’elle continue a fonctionner avec des fonctions admettant une singularite integrable en une borne d’integration. Supposons par exemple que f est de la forme
sur I’intervalle [0, I], avec /L de classe C”, P(X) =
J; t2”-2(1
Si
t27n-2(1
h(O) # 0 et 0 < n: < 1. En posant - t)27rL-2dt avec ,~12 > ~ - > t)2nL-2&
_
(11)
on a g(u) = f(P(U))P’(,U)
= h(pjJg’(U).
par consequent, la dCrivCe d’ordre 2k de g(u) est Cquivalente en 0 a XU~‘~‘-~-~(~~‘-‘)-‘~. Au vue de la formule (4), et a condition que 771 soit assezgrand, on pourra obtenir une methode d’ordre 2k. Cette idee est precide dans le theoreme suivant TH~OR~ME 2. - Soit f une application comme dans (10) et P le polyn8me dkjini par (1 1 ), alors il existe une constante C, > 0 telle que
IW)
- ~df)l
5 GN(2m~l~(l+-F
pour tout E > 0.
(12)
Si la singularite’ en 0 est de type logarithmique, alors il suffit de remplacer IZIpar 0 dans l’estimation ci-dessus. Demonstration. - On se ramene d’abord par un developpement limit6 en 0 a l’integration de l/gF. D’autre part on a: 1 - cos(27rZNt) = 2sin2(rZNt) - 2(z-lNt)2 au voisinage de 0, et la fonction f(2”)(t) 1-cos(27rZNt) sera integrable sur [O:l] si et seulementsi l’exposant de t au voisinage de 0, qui (27rZN)2” est ici 2m - 2- o(2rn - 1) - 2k + 2 est > -1. Si la partie entiere de 2m( 1 - (Y) + a + 1 est paire, on peut done poser 2k = E( (2 m - l)( 1 - a)) + 2 oti E(-) designe la fonction partie entiere. Posonsalors 2j? = 2F - (277~- l)(l - a) + E. On peut choisir E > 0 assezpetit tel que 0 < 2,6 < 2. On a alors
I
fc2”‘@)
846
I
1 - cos(27rlNt) (zTzN)2k
M&ration
numkrique
de fonctions
singulikes
ou
non
pour E > 0: ]f(‘“)(t)l]t]‘” est dans L’(]O, l[) ce qui prouve l’estimation. Si la partie entiere de 2n~(l - o) + cy + 1 est impaire, un raisonnement analogue conduit au meme resultat. Pour la singularite logarithmique, la demonstration est identique. 0 4. Applications 4.1. Ordre
numkriques
de la mkthode
de pkriodisation
pour des fonctions
reguli&res
Dans ce paragraphe, now appliquons la mtthode de periodisation pour diverses valeurs de k et pour la fonction f(z) = exp(r) sur I’intervalle [0, 11. Nous avons fait varier le nombre de points d’integration N de 100 a 300, notre but Ctant de verifier que l’estimation d’erreur (9) est optimale. Nous cherchons pour chaque A; a montrer numeriquement que
et a determiner Ck et y. Cette determination se fait en traqant le logarithme de I’erreur en fonction de log N et en Cvaluant la pente de la courbe obtenue. Pour &tre exploitables, ces calculs ont dQ &tre realids avec une precision de plusieurs centaines de chiffres significatifs. On obtient les resultats suivants
1
C’k
k
2k
II(f) - fflO(.f)l
V(f) - ffm(.f)I
2
1.S571
x 10”
3.999s
4
O.lS39.5
x lo-:’
3
3.7144
x 10’
5.999s
c
0.36831
x lo-
4
1.1098
x 10:’
7.99ss
8
0.10514
x lo-’
0.42957
x 10-i
5
4.S.341
x 10-k
9.99Gl
10
0.48090
x lo-’
0.45S2.5
x lo-’
10
3.6130
x 10”
19.976
20
O.G27OG
x lo-:$
o.s1;9s
x lo-”
’
0.113SS
X lo-’
0.57934
x lo-”
L’estimation d’erreur que nous avons obtenue est done clairement optimale. On observe une croissance importante des constantes CI, quand k augmente. Cette mtthode a pour vocation d’etre utilisee sous sa forme composite et done en general, N sera compris entre 5 et 20. Pour ces valeurs de N, un bon choix de k semble ttre 5 1. k 5 10. 4.2. Ordre de la mkhode de pkriodisation pour des fonctions admettant une singularit
en 0
Nous effectuons les mtmes calculs qu’au paragrapheprecedent, mais cette fois-ci pour fi (T) = & et f2(2) = in(z). Nous obtenons les resultats suivants pour la fonction fl
U1
cm
i
(2ru
- l)(l
- 0)
V(f)
-
ffl(l(f)l
V(f) - ff%fl(f)l
2
1.215s
2.0101
2
O.llG51
x 10-I
0.29033
x lo-
3
0.85383
3.3333
3.3333
0.31289
x lo-”
0.3341S
x lo-”
3
39.567
6.0606
G
0.24947
x lo-”
0.4.5152
X 1OF
O.GS427
x lo-”
0.13577
x lo-!’
10
0.46417
x 10”
12.400
12.GG7
847
P. Helluy
et al.
et les rksultats
suivants
pour la fonction
f2
(21,) - 1)
P(f)
- HlO(.f)l
I??
cm
2
3.2220
3.0303
3
0.30240
x 10-z
3
34.022
5.0503
5
0.30379
X
5 10
30485.0
9.0909 19.06i
9 19
0.14231 04i492
1.734T
x 10”’
L’estimation d’erreur ( 12) semble done optimale. Pour des valeurs de N de I’ordre de 10, un bon choix de
m
IO-.’
II(f)
- fboC.f)l
0.36694
x lo-:’
0.033x
x 1o-5
x lo-’
0.42332
x lo-’
x lo-’
0.91853
X ltl-”
est ici aussi 5 5 n-r, 5 10.
Refkrences bibliograpbiques [l]
Davis P., Rabinowitz P.. Methods of numerical integration, second edition, Comp. Sci. and Applied Math., Academic Press, 1984. [2] Helluy P., Maire S., Ravel P., (to appear). 131 hi M., Moriguti S., Takazawa Y., On a certain quadrature formula (in Japanese), Kokyuroku of Res., Inst. for Math. Sci., Kyoto Univ., Vol. 91, 1970, pp. 82-l 18. [4] Murota K., Iri M., Parameter tuning and repeated application of the IMT-type transformation in numerical quadrature, Numer. Math. 38 (1982) 347-363. [5] Takahasi H., Mori M.. Quadrature formulas obtained by variable transformation, Numer. Math. 21 (1973) 206219.
848