Solução Numérica para Equações Integrais com Núcleo
Fracamente Singular.
Leonardo Guillermo Felipe1
1
Centro de Engenharias e Ciências Exatas da Universidade Estadual do Oeste do
Paraná.
Rua da Faculdade, 645- Jd. Santa Maria-Toledo-PR
[email protected]
Resumo: O objetivo deste trabalho é desenvolver e discutir a convergência do
método numérico integração produto para aproximar a solução das equações
integrais de Volterra de segunda espécie com núcleo fracamente singular. Testes
numéricos são realizados, os quais confirmam o resultado teórico sobre a taxa de
convergência do método numérico proposto.
Palavras chaves. Convergência. Integração produto. Quadratura.
1. Introdução.
Consideremos a equação integral de Volterra de segunda espécie,
x
y ( x) = g ( x) + ò p( x, s ) K ( x, s , y ( s )) ds ,
a £ x < ¥,
(1)
p
g
onde o núcleo
é fracamente singular e as funções
e K são dadas e assumidas
suficientemente suaves de modo a garantir a existência e unicidade da solução y Î C[ a
,b], ver por exemplo, [ATKINSON, 1974] e [BRUNNER , HOUWEN,1986].
Formas típicas para p (x,s) são,
a
p( x, s ) = x - s
-l
, 0 < l <1
.
Equações do tipo (1) surgem no estudo de problemas de condução de calor,
superfluidez e eletroquímica (ver por exemplo, [KELLER , OLMSTEAD,1971-2] ,
[LEVINSON, 1960].
Para equações integrais de Volterra com núcleo limitado, a suavidade do núcleo e
da função g (x) determinam a suavidade da solução no intervalo [ a ,X] para X> a . Por
outro lado, se considerarmos o núcleo como sendo fracamente singular, então a solução
resultante é tipicamente não suave no ponto inicial do intervalo de integração, onde sua
derivada torna-se não limitada. Alguns resultados relacionados com o comportamento
da solução exata das equações integrais do tipo (1) podem ser encontrados em
[BRUNNER , HOUWEN, 1986], e [LUBICH, 1983].
Métodos numéricos para resolver equações integrais de tipo Abel de segunda
espécie, isto é, do tipo (1) tem sido estudados por [BRUNNER , HOUWEN, 1986],
[CAMERON , McKEE,1984], [RIELE, 1982], [ABDALKHANI, 1990] e [LUBICH,
1985], entre outros autores.
O método para a solução da equação integral (1) que consideramos neste trabalho
procura superar a dificuldade causada pelo mau comportamento da solução y (s) no
ponto inicial s= a . Assim, dado um intervalo relativamente pequeno [ a ,b], primeiro
resolvemos o problema,
x
y ( x) = g ( x ) + ò p ( x, s ) K ( x, s, y ( s )) ds,
a £ x £ b,
(2)
pelo método de Nystrom sobre todo o intervalo, o qual é baseado na regra de integração
produto de tipo interpolatório. Após o intervalo inicial, o mau comportamento da
derivada de y é de menor importância. Logo, resolvemos o problema,
a
x
y ( x) = g1 ( x ) + ò p ( x, s ) K ( x, s, y ( s )) ds,
b
b£x<¥
,
(3)
onde,
b
g1 ( x) = g ( x ) + ò p ( x, s ) K ( x, s, y (s )) ds ,
a
por um método padrão passo a passo para solução regular. Uma vez que os cálculos de
g1(x) dependem das aproximações iniciais de y(x) para x Î [ a ,b], os dois métodos
devem ser considerados como complementares.
Na seção 2 deste trabalho descrevemos a parte inicial do método para a equação
integral (2) e na seção 3 serão dados resultados de convergência uniforme. A aplicação
do método inicial junto com o esquema padrão para a solução da equação integral (3)
serão fornecidos na seção 4. Nossos resultados numéricos serão comparados com os
resultados disponíveis na literatura obtidos por métodos alternativos.
2. O Método de Integração Produto.
Nesta seção descrevemos o método de Nystrom para resolver numericamente a equação
(2). Por conveniência e sem perda de generalidade, assumiremos que [ a ,b]=[-1,1],
x
y ( x) = g ( x ) + ò p ( x, s ) K ( x, s, y ( s )) ds,
-1
- 1 £ x £ 1.
(4)
Tendo escolhido N+1 pontos distintos {xn}, para n=0,...,N, no intervalo [-1,1],
colocamos a equação (4) nos nós {xn} :
xn
y ( x n ) = g ( x n ) + ò p( x n , s ) K ( x n , s , y ( s )) ds, n = 0,1,....., N .
-1
Logo, usamos os polinômios de interpolação de Lagrange,
N
L N ( K ; s ) = å l N , j ( s ) K ( x n , x j , y ( x j ))
j =0
para aproximar K (xn,s,y(s)) e obter o seguinte:
N
y N , n = g ( x n ) + å wn , j K ( x n , x j , y N , j ),
j =0
onde,
n = 0,1,......, N ,
(5)
xn
wn, j = ò p ( x n , s ) l N , j ( s ) ds
-1
.
3. O Método Inicial: Convergência.
Em toda esta seção, o símbolo C denotará uma constante positiva. Para fazer a análise
de convergência, usaremos a equação teste seguinte,
x
y ( x) = g ( x ) + ò p( x, s ) y ( s ) ds ,
- 1 £ x £ 1,
-1
(6)
p
onde assumiremos que a função g Î C[-1,1] e o núcleo
é fracamente singular.
Consequentemente, a equação (6) tem uma única solução y Î C[-1,1] com,
possivelmente, derivada não limitada no ponto x= -1.
Considerando uma malha {xj} para j=0,...,N e aplicando o método dado em (5)
para a equação teste (6), a solução aproximada y N(x) será obtida mediante o método de
Nystrom,
N
y N ( x) = g ( x) + å w j ( p, x) y N ( x j ),
j =0
onde,
x
w j ( p, x) = ò p( x, s ) lN , j ( s ) ds
-1
.
O método numérico é convergente de ordem r em [-1,1] se e somente se, para N
suficientemente grande existe uma constante C tal que,
y ( x ) - y N ( x) ¥ £ CN - r
.
Para fazer a análise de convergência uniforme da solução aproximada yN(x) para
solução exata y(x) da equação (6), observamos que,
y ( x) - y N ( x) = å w j ( p, x ) { y ( x j ) - y N ( x j )} + t N ( p, y , x)
N
j=0
onde tN( p ,y,x) é o erro de truncamento local definido por,
n
t N ( p, y, x ) = ò p ( x, s ) y ( s ) - å w j ( p, x) y ( x j ).
x
-1
j =0
Logo,
y - yN
¥
£ ( I - AN ) -1
¥
tN
¥
,
(7)
onde AN é o operador linear definido por,
N
AN f ( x) = å w j ( p, x) f ( x j ), f Î C[-1,1], x Î [-1,1]
j =0
(8)
Antes de prosseguir, estabelecemos um resultado sobre as propriedades de convergência
da regra de quadratura-produto:
Teorema 1. Seja {xj}, j=0,...,N, os zeros de um conjunto de polinômios de grau N+1 que
são ortogonais sobre [-1,1] em relação a função peso,
w( s ) = u ( s ) (1 - s )a (1 + s ) b , - 1 < a £ 3 / 2, b ³ -1 / 2.
(9)
onde, u (x) é positiva e contínua em. Seja LN(f;s) o polinômio de interpolação de grau £
N que coincide com a função f nos nós {xj}. Então, para toda função f contendo
s
singularidade do tipo (1+s) , s >-1 (não inteiro), e em particular para toda função fÎ
C[-1,1] é verdade que,
lim
N ®¥
ò
x
-1
x
p ( x, s ) f ( s ) ds - ò p( x, s ) LN ( f , s )ds
-1
=0
¥
(10)
Em particular, temos o limitante,
t N ( x - s ) - l , f , x = O{( N + 1) -2 - 2s + 2l log( N + 1), 0 < l < 1, s > -1}
¥
( I - AN )-1
Agora analisaremos o comportamento de
propósito, estabeleceremos alguns lemas preliminares.
¥
(11)
em (7). Para este
Lema 1. Seja o conjunto de nós {xj} como definido no teorema 1 com a restrição -1/2<
a , b <3/2. Seja l (s) o j-ésimo polinômio fundamental de Lagrange . Para [c,d] Í [N,j
1,1] existe um número C positivo e q>1 tal que,
N
sup å
N
j =0
ò
d
c
d
q
p (d , s ) l N , j ( s ) ds £ C é ò p (d , s ) ds ù
êë c
ûú
1/ q
para todo p Î Lq com,
p
q
q
é 1
ù
= ê ò p (d , s ) ds ú
-1
ë
û
(12)
1/ q
.
Lema 2. Com as mesmas hipóteses do lema 1 e com o núcleo p satisfazendo,
p Î Lq , q > 1,
ì
ï
í
lim p( x ' , s ) - p( x, s ) = 0
ï
q
îx ' ® x
para todo x Î [-1,1] . Então,
(13)
N
lim
sup å w j ( p , , x , ) - w j ( p, x) = 0
,
x ®x
N
j =0
(14)
para todo x Î [-1,1] , onde p´= p(x´,s).
Teorema 2. Seja AN o operador definido em (8) e os nós {xj} como no Lema 1. Se (10),
(12) e (14) são verdadeiros, então para todo N suficientemente grande existe uma
constante C>0 independente de N tal que,
( I - AN ) -1 £ C
¥
.
Lema 3. O núcleo p satisfaz a condição (13).
Teorema 3. Seja y a solução exata da equação (6). Seja yN a solução aproximada obtida
pela discretização do termo integral da equação (6) usando a regra de quadraturaproduto do tipo interpolatório construída sobre um conjunto de nós {xj}. Se os nós {xj}
os zeros de um conjunto de polinômios de grau N+1 que são ortogonais sobre [-1,1] em
relação a função peso (9) com -1/2< a , b <3/2; então, y converge uniformemente para
N
y. Além disso, a taxa de convergência de yN para y coincide com a da regra de
quadratura escolhida para aproximar o termo integral na equação (6).
4. Resultados Numéricos e Discussâo.
O método (5) baseado nos nós de Radau ( isto é, os nós coincidem com os zeros dos
P (1,0 ) ( x)
) usando a regra de Simpson-produto sobre malha
polinômios de Jacobi N
uniforme já foi implementado em [CAMERON , McKEE, 1984] para resolver as
equações seguintes (tomadas de uma coleção de problemas propostos em [BRUNNER ,
HOUWEN, 1986] e [RIELE, 1982],
1 x
y ( x) = 1 - e - x ( x - s ) -1 / 2 y (s ) ds
ò
0
p
,
(15)
1 x
y ( x) = e erfc ( x ) - e - x + 2 p F ( x )
2
,
{
onde,
F (t ) = e -t
}
2
t
òe
0
u2
du
;
1
p 1
æ1- x ö 1 x
-1 / 2
+ - arcsenç
y ( s ) ds
÷ - ò0 ( x - s )
8
4
1
+
x
4
1+ x
è
ø
,
1
y ( x) =
1+ x .
y ( x) =
(16)
Para h dado, fixamos o ponto b no intervalo ( a ,X), com X> a . Logo, resolvemos
o problema (2) pelo método (5) baseado em N+1 nós de Radau no intervalo [ a ,b]. Se a
s
solução y(s) contém uma singularidade da forma (s- a ) , s >-1; então, a ordem de
-2 - 2s + 2 l
log( N + 1) conforme foi estabelecido em (11)
convergência do método é ( N + 1)
do teorema 1. Com b fixado no intervalo [b,X] definimos a malha x o º b,
xn = x0 + nh, n = 1,2,..., M ; xM º X e logo aplicamos o método de Simpson-produto, ver
[CAMERON , McKEE, 1984]. Os valores aproximados para yn, n=2,...,M são obtido de,
ìy ( x ) = g ( x ) + xn ( x - s ) -l K ( x , s , y ( s )) ds,
1
n
n
òx0 n
ï n
í
x
ïg1 ( x n ) = g ( n ) + ò 0 ( x n - s ) -l K ( x n , s , y ( s )) ds, l £ s £ x 0
a
î
.
y
Desde que a solução (x) é suave no intervalo [b,X], (ver [LUBICH, 1983] ), este
4-l
método converge do mesmo modo que h ( ver [CAMERON , McKEE, 1984] ). Os
valores iniciais y0 e y1 são determinados usando o método de Nystrom no intervalo [ a
,x0] e a regra passo a passo de Simpson-produto é aplicada nos pontos x0-h, x0, x0+h.
Para testar o desempenho do método apresentado neste trabalho, comparamos os
nossos resultados numéricos obtidos com os correspondentes resultados existentes na
literatura. Assim, damos especial atenção aos resultados fornecidos por [HAIRER,
LUBICH , SCHLICHTE, 1988] que utilizaram o método de passo múltiplo introduzido
por [LUBICH, 1985], o qual é considerado um método eficiente para equações integrais
de tipo Abel, segundo as considerações de [BAKER, 1987].
Na Tabela 1 denotamos por IP os resultados numéricos deste trabalho, por P o método
inicial (5), por L os resultados numéricos de [HAIRER, LUBICH , SCHLICHTE, 1988]
e por BE os resultados em relação a equação (15) apresentados por [BRUNNER ,
EVANS, 1977], que utilizaram o método de colocação com polinômios por pedaços e
malha uniforme. O método IP foi aplicado com b=0,2 e h=0,01.
Tabela 1. Erros absolutos para a solução aproximada da equação (15)
em x=1 calculados por métodos diferentes
N
4
5
8
10
16
20
32
IP
8,7E-8
1,0E-8
2,1E-9
8,7E-10
1,0E-10
3,6E-11
3,8E-12
P
4,3E-5
1,8E-5
7,3E-7
1,2E-7
6,6E-9
8,6E-10
1,1E-10
BE
L
3,4E-8
8,4E-4
4,1E-8
2,3E-8
5,0E-5
4,7E-9
Em relação a equação (16), a Tabela 2 mostra os resultados fornecidos por :
1. [LINZ, 1969], denotado por LZ, utilizou o método de integração produto de ordem
p=3 baseado na regra de Simpson generalizado.
2. [EL TOM, 1974], denotado por ET, utilizou funções spline de grau 5.
3. [CAMERON , McKEE, 1984], denotado por CK, utilizou um método de integração
produto de tipo interpolatório.
Salientamos que os resultados de convergência fornecidos pelos métodos de LZ, ET e
CK são válidos se a solução y(x) é suficientemente suave.
O método deste trabalho (IP) foi implementado com b=0,02 e h=0,002.
Tabela 2. Erros absolutos para a solução aproximada da equação (16) em x=1
calculados por métodos diferentes.
N
4
5
8
10
16
20
32
64
100
128
200
256
IP
1,4E-9
1,9E-10
9,7E-13
P
2,8E-6
3,1E-7
6,3E-10
1,2E-11
1,3E-16
LZ
ET
£ 7,0E-8
CK
L
3,0E-7
5,6E-6
£ 2,0E-9
1,3E-7
8,0E-7
1,4E-8
1,0E-7
1,9E-9
5,6E-10
2,8E-14
6,7E-11
5,6E-16
5,7E-12
Os resultados numéricos fornecidos acima mostram que o método apresentado
neste trabalho confirma a taxa de convergência teórica.
No caso da equação (15) (ver Tabela 1) o método IP, com [ a ,b]=[0;0,2], h=0,01. e
o método L envolvem N=8 e N=32 nós, respectivamente, para encontrar a solução
aproximada em x=1 com erro absoluto na ordem de 10-9.
O comportamento dos erros na Tabela 2 é devido ao fato que a solução da equação
(16) é suave no intervalo de integração. Em particular, podemos observar que, se a
equação (16) for resolvida em 0 £ x £ 0,02 pelo método inicial P; o método IP, com
h=0,002 requer N=8 para encontrar a solução aproximada em x=1, com erro absoluto na
ordem de 10-12 . A mesma exatidão é obtida pelo método L em x=1 com N=256.
5. Conclusões.
Neste trabalho foi analisado a convergência do método integração produto para
equações integrais de Volterra com núcleo fracamente singular. Assim, se considerarmos
o núcleo fracamente singular, então a solução resultante é tipicamente não suave no
ponto inicial do intervalo de integração, onde sua derivada torna-se não limitada. Para
amenizar esta dificuldade, temos proposto o método de integração produto junto a um
método passo a passo apropriado para a solução regular. Conclui-se que esta técnica
pode ser uma ferramenta numérica importante para aproximar a solução de equações
integrais de Volterra de segunda espécie com núcleo fracamente singular. Os testes
numéricos realizados confirmam o resultado teórico de convergência fornecida neste
trabalho.
6. Referências.
ABDALKHANI, J. A numerical approach to the solution of Abel integral equations of
the second kind with nonsmooth solution. J. Comput. Appl. Math., v. 29, p. 249-255.
1990.
ATKINSON, K.E. An existence theorem for Abel integral equations. Siam J. Math.
Anal., v.5, p.729-736. 1974.
BAKER, C.T.H. The state of the art in the numerical treatment of integral equations,
The state of the Art in Numerical Analysis. Clarendon Press, Oxford. 1987.
BRUNNER, H.; EVANS, M.D. Piecewise polynomial collocation for Volterra-Type
integral equations of the second kind. J. Inst. Math. Appl. v. 20, p. 415-423. 1977.
BRUNNER, H. ; Te RIELE, J.J. Volterra-type integral equations of the second kindwith
nonsmooth solutions. J. Integral Equations, v.6, p.187-203. 1984.
BRUNNER, H. ; HOUWEN, P.J. The numerical solution of Volterra equations. NorthHolland, Amsterdam. 1986.
CAMERON, R.F. ; McKEE, S. Product integration methods for second kind Abel
integral equations. J. Comput. Appl. Math., v.11, p.1-10. 1984.
EL TOM, M.E.A. Spline function approximations to the solution of singular Volterra
integral equations of the second kind . J. Inst. Math. Appl. v.14, p.303-309. 1974.
HAIRER, E.; LUBICH, Ch. ; SCHILCHE, M. Fast numerical solution of weakly
singular Volterra integral equations. J. Comput. Appl. Math., v.23, p.87-98. 1988.
KELLER, J.B. ; OLMSTEAD,W.E. Temperature of a nonlinear radiating semiinfinite
solid.Q. Appl. Math. v.29, p.559-566. 1971-72.
LEVINSON, J.E. A nonlinear Volterra equation arising in the theory of superfluidity. J.
Math.Analysis Appl. v.1, p.1-11. 1960.
LINZ, P. Numerical Methods for Volterra integral equations with singular kernels. Siam
J. Numer. Anal. v.6, p. 365-374. 1969.
LUBICH, CH. Runge-Kutta theory for Volterra and Abel integral equations of the
second kind. Math. Comp. v.41, p. 87-102. 1983.
RIELE, H.J.J. Collocation methods for weakly singular second kind Volterra integral
equations with non-smooth solution. IMA J. Numer. Anal. v.2, p. 437-449. 1982.