Embed Size (px)
Citation preview
Volume 1 – N. 1 – 2012
Solução Numérica de Equações Diferenciais de Ordem Fracionária aplicadas em
Modelo de Compartimentos Maysa Costa de Castro, Fernando Luiz Pio dos Santos ...................................................02
Os Códigos Controle da Paridade CPq(n) obtidos por restrição de um Código de
Goppa Racional Jaime Edmundo Apaza Rodriguez ...............................................................................09
Cálculo Fracionário Aplicado ao Problema da Tautócrona Pedro Felipe Pavanelo Ramos, Rubens de Figueiredo Camargo .......................................15
Métodos Iterativos para a Solução da Equação de Poisson Valdirene da Rosa Rocho, Dagoberto Adriano Rizzotto Justo...........................................23
Modelo dinâmico incluindo Manejo Integrado de Pragas (MIP) e estrutura espacial
no combate àDiaphorina Citri Priscila Azevedo da Silveira .......................................................................................28
Solucao Numerica de Equacoes Diferenciais de OrdemFracionaria aplicadas em Modelo de Compartimentos
Maysa Costa de Castro, Fernando Luiz Pio dos Santos ∗
31 de outubro de 2012
Resumo
O proposito deste trabalho foi o estudo de modelagem matematica e resolucaonumerica utilizando modelos de compartimentos. Realizou-se medidas de concen-tracao de substancias em um dado sistema fısico fazendo uso de equacoes diferenciaisde ordem fracionaria. Os resultados numericos foram obtidos atraves da imple-mentacao das equacoes modelo considerando um sistema de dois compartimentos.Uma comparacao entre as solucoes numerica e analıtica foi efetuada. Evidenciou-seo decaimento da concentracao da substancia ao longo do tempo em ambos compar-timentos.
Palavras Chave: Equacoes diferenciais de ordem fracionaria, distribuicao da con-centracao, solucao analıtica.
1 Introducao
Na literatura pode-se encontrar varias definicoes alternativas sobre derivadas de or-dem fracional. Conforme em [9], serao consideradas as definicoes de Caputo e deGrunwald-Letnikov. Pela definicao de Caputo para derivadas de ordem fracionarias,tem-se:
aDαt f(t) =
1
Γ(n− α)
t∫a
fn(τ)
(t− τ)(α− n+ 1)dτ (1.0.1)
sendo n − 1 < α < n a ordem fracionaria. Pela definicao de Grunwald-Letnikov,tem-se que:
aDαt f(t) = lim
h→0
1
hα
[k]∑j=0
(−1)j(
αj
)f(t− jh) (1.0.2)
em que [k] e a parte inteira de k = t−ah , com a e t sendo os limites reais do oper-
ador aDαt ; h e o passo-no-tempo; A seguinte expressao para o calculo do coeficiente
∗Email: [email protected] , [email protected]
02
binomial cαj = (−1)j(
αj
), j = 0, 1, ..., pode ser utilizada:
cαj = (1− 1 + α
j)c
(α)j−1, (1.0.3)
assumindo cα0 = 1. A solucao geral numerica da equacao diferencial fracionaria
aDαt = f(y(t), t) e dada por:
y(tk) = f(y(tk), tk)hα −
k∑j=1
cαj y(tk−j), k = 1, ... (1.0.4)
Na proxima secao e descrita a modelagem para a obtencao dos resultados numericosda distribuicao da concentracao em um sistema de compartimentos.
2 Metodologia
O esquema do modelo de dois compartimentos [4, 9, 7, 10] utilizado neste trabalhoe ilustrado na Figura (1).
Figura 1: Modelo de dois compartimentos.
Pelo metodo de solucao geral descrito pela Equacao (1.0.4), obtem-se o seguintemodelo numerico para o sistema de compartimentos anterior, Figura (1):
q1(tk) = (−K21q1(tk−1))hα1 +
k∑j=v
cα1j q1(tk−j) (2.0.5)
q2(tk) = (K21q1(tk)−K02q2(tk−1))hα2 +
k∑j=v
cα2j q2(tk−j) (2.0.6)
onde K21 e K02 sao as taxas de transferencias do compartimento 1 para o 2 e de 2para o exterior, respectivamente; α1 = α2 = α e a ordem fracionaria.
O algoritmo para o calculo da concentracao de uma substancia ao longo dotempo em 1 e 2 segue os passos:
• Passo 1: Entrada dos parametros α, K21, K02 e d mg/L; Prescricao dascondicoes iniciais d1 = 100 mg/L (dose inicial), K01 = K12 = 0, q1(t0) = d1 eq2(t0) = 0;
• Passo 2: Calculo dos coeficientes binomiais pela Equacao (1.0.3);
03
• Passo 3: Calculos de q1(tk) e q2(tk) utilizando as Equacoes (2.0.5) e (2.0.6);Finalmente,
• Passo 4: Obtecao das concentracoes C1 e C2:
C1 = q1/v (2.0.7)
C2 = q2/v (2.0.8)
sendo v = d1/d o volume dos compartimentos.
3 Resultados
Os resultados numericos foram obtidos implementando-se o algoritmo descrito nasecao anterior no programa MatLab. A Tabela 1 mostra os parametros utilizadospara os testes numericos.
Tabela 1: Parametros.Teste α K21 K02 d(mg/L)1 0.93 2.89 2.74 30.172 0.60 0.90 1.68 35.003 1.00 0.34 0.37 4.00
Os graficos abaixo mostram as concentracoes nos compartimentos utilizando osparametros da Tabela 1
Figura 2: Teste 1: Decaimento da concentracao nos compatimentos 1 e 2.
O Teste 1 evidencia o decaimento da concentracao da substancia em ambosos compartimentos ao longo do tempo, conforme pode ser visto na Figura (2).Essa substancia e transferida de 1 para 2 com taxa K21. Observa-se nıveis deconcentracao distintos em cada compartimento, sendo mais elevada em 1, onde adose inicial injetada. Alem disso, e possıvel notar o comportamento de decaimentoexponencial. A analise no compartimento 2 dos demais Testes, 2 e 3, tambem mostraeste aspecto da solucao, Figuras (3) e (4).
04
0 5 10 15 20 250
0.2
0.4
0.6
0.8
1
1.2
1.4
Tempo (h)
Con
cent
rção
(m
g/L)
Compartimento 2
Figura 3: Teste 2: Decaimento da concentracao no compartimento 2; α = 0.60.
0 5 10 15 20 250
0.5
1
1.5
Tempo (h)
Con
cent
rção
(m
g/L)
Compartimento 2
Figura 4: Teste 3: Decaimento da concentracao no compartimento 2; α = 1.00.
Com o intuito de se verificar o correto comportamento da solucao encontradapelo modelo de equacao de ordem fracionaria, comparou-se com a solucao analıticaobtida modelando-se o sistema de compartimento pela seguinte equacao matricial:
C′(t)−KC(t) = 0 (3.0.9)
sendo C(t) a concentracao no tempo t; K a taxa de transferencia. Com a solucaoapresentada abaixo:
C(t) = C0eKt (3.0.10)
sendo C(t) =
(C1
C2
), C0 =
(C0
)e K =
(−k21 0k21 −k20
). Com −k21 e
−k02 sao os autovalores de K.
Com isso podemos dizer que
K = H
(e−k21t 0
0 e−k02t
)H−1 (3.0.11)
ondeH =
(−k21 − k02 0
−k21 −k21
)e sua inversıvel e dada porH−1 =
( 1k21−k02
0
− 1k21−k02
− 1k21
)
05
Obtendo-se assim a solucao analıtica geral, dada pelo sistema de equacoes abaixo:C1 = C0e
−K21t
C2 =K21C0
K21−K02[e−K02t − e−K21t]
sendo C0 a dose, K21 a taxa de transferencia de 1 para 2; K02 a taxa de transferenciade 2 para o exterior.
A analise do sistema foi realizado utilizando os parametros do Teste 3 (ordem fra-cionaria α = 1), mostrando que as curvas da solucao analıtica, dada pela Equacao(3), e da solucao numerica, obtida por (2.0.7), encontram-se bastante proximas,como mostra a Figura (5) (a). Esta proximidade entre as solucoes analıtica enumerica pode ser verificada pelo refinamento do intervalo das solucoes, conformeilustra a Figura (5) (b). O erro relativo entre as solucoes para um subintervalo detempo foi calculado e mostrou-se pequeno. Os resultados deste calculo podem servisualizados na Tabela 2 para o compartimento 1 e na Tabela 3 para o comparti-mento 2.
0 5 10 15 20 250
0.5
1
1.5
2
2.5
3
3.5
4
Tempo (h)
Con
cent
raçã
o (m
g/L)
Solução analítica do compartimento 1Solução numérica do compartimento 1Solução analítica do compartimento 2Solução numérica do compartimento 2
(a)
3.065 3.07 3.075 3.08 3.085 3.09 3.095 3.1 3.105
1.401
1.402
1.403
1.404
1.405
1.406
1.407
1.408
Tempo (h)
Con
cent
rção
(m
g/L)
Solução analítica do compartimento 1Solução numérica do compartimento 1Solução analítica do compartimento 2Solução numérica do compartimento 2
(b)
Figura 5: Teste 3: (a) Comparacao e (b) Erro entre as solucao numerica e analıtica para ambos
os compartimentos; α = 1.00.
06
Tabela 2: Teste 3: Comparacao entre as solucoes numerica e analıtica da concentracaono compartimento 1.
Tempo Analıtica Numerica Erro relativo(%)3.073 1.407 1.408 0.073.075 1.406 1.407 0.073.083 1.402 1.403 0.073,085 1.401 1.402 0.07
Tabela 3: Comparacao entre as solucoes numerica e analıtica da concentracao no com-partimento 2.
Tempo Analıtica Numerica Erro relativo3.063 1.405 1.404 0.073.083 1.404 1.403 0.073.105 1.403 1.402 0.07
4 Conclusao
Nos testes realizados observou-se nıveis de concentracao distintos em cada compar-timento, mais elevada no primeiro, onde a dose inicial e injetada. Alem disso, odecaimento das concentracoes no intervalo apresentam um comportamento expo-nencial. A solucao analıtica para o caso em que a ordem da equacao e inteira foicalculada e comparada com a respectiva solucao numerica. O intuito foi de verificaro correto comportamento do decaimento exponecial da concentracao nos comparti-mentos. As curvas mostraram-se bem proximas. O erro relativo entre as solucoespara este caso foi calculado e mostrou-se pequeno. Concluı-se, portanto, que asequacoes de ordem fracionaria aplicadas no modelo de dois compartimentos foi ca-paz de evidenciar corretamente o compartamento da distribuicao da concentracaoao longo do tempo.
Referencias
[1] BAILEY M. J., SHAFER L. S., A simple analytical solution to the three-compartment pharmacokinetic model suitable fo computer-controlled infusionpumps, IEEE Transactions on biomedical engineering 38,(1991) 522-525.
[2] BOYCE, W.E. - DIPRIMA R.C.- Elementary differential equations, JohnWiley, New York, 1969.
[3] CAMARGO F. R.,- Calculo fracionario e aplicacoes. Tese de doutorado.Universidade Estadual de Campinas, Campinas(2009).
[4] COBELLLI C., FREZZA M., TIRIBELLI C., Modeling Identification andParameter Estimation of Bilirubin Kinetics in normal, Hemolytic and Gilbert’s
07
States. Computers and Biomedical Research 8,(1975) 522-537.
[5] FIGUEIREDO D. G.- Equacoes diferenciais aplicadas (3edicao), Rio deJaneiro: IMPA, 2007.
[6] HOLMES H. M., Introduction to NUmerical Methods in DifferentialEquations- Springer: 2007.
[7] KIBBY R. M., The matrix representation pharmacokinetic models, J. theorBiol 77,(1979) 333-348.
[8] LAWRENCE X. Y., GORDON L. A., A compartmental absorption andtransit model for estimating oral drug absorption. Elsevier Internation Journalof Pharmaceutics 186,(1999) 1119-125.
[9] PETRAS I., MAGIN R. L., Simulation of drugs uptake in a two compartmen-tal fractional model for a biological system.Elsevier Commun Nonlinear SciNumer Simulat 16,(2011) 4588-4595.
[10] SAVIC M. R., JONKER M. D., KERBUSCH T., KARLSSON O. M.,- Im-plementation of a transit compartment model for describing drug absorptionin pharmacokinetic studies, J. Pharmacokinet Phamarcodyn (2007)34:711-716.
[11] UPTON N. R., The two-compartment recirculatory pharmacokinetic model -an introduction to recirculatory pharmacokinetic concepts, British Journal ofAnaesthesia 92(4),(2004) 475-484.
08
Os Codigos Controle da Paridade CPq(n) obtidos porrestricao de um Codigo de Goppa Racional
Jaime Edmundo Apaza Rodriguez ∗
Departamento de Matematica, UNESP, Ilha Solteira
18 de outubro de 2012
Resumo
Os Codigos Controle da Paridade surgiram como um exemplo de uma famılia decodigos detectores de um erro simples. Seu nome provem de um codigo binario queacrescentava um sımbolo extra para que o numero de 1’s fosse par. Em 1977, V.D. Goppa introduziu uma nova forma de construir codigos lineares usando curvasalgebricas definidas sobre corpos finitos. Esses codigos sao chamados de CodigosGeometricos de Goppa. Neste trabalho usamos algumas das ideias de Goppa paraconstruir uma certa classe de Codigos de Controle da Paridade CPq(n) sobre o corpoFq. Construimos o Codigo de Controle da Paridade CPq(n) como restricao de umCodigo de Goppa Racional. No processo de construcao sao considerados dois casos,dependendo de se o comprimento n do codigo e divisıvel ou nao pela caracterısticado corpo finito Fq.
Palavras Chave: Codigos Lineares, Corpo Finito, Curva Algebrica nao-singular,Codigos Geometricos de Goppa, Teorema de Riemann-Roch.
1 Introducao
Os Codigos Controle da Paridade surgiram no seculo passado como exemplo de umafamılia de codigos detectores de um erro simples. Seu nome provem de um codigobinario que acrescentava um sımbolo extra para que o numero de 1’s fosse par.
Os Codigos Controle da Paridade sao utilizados pelo fato de usar um numeromınimo de sımbolos de redundancia, o que reduz custos de implementacao e deco-dificacao em sistemas de erro simples.
Os codigos de Goppa ilustram uma construcao pratica para famılias de bonscodigos que ultrapassam determinadas cotas assintoticas conhecidas, mostrando seruma boa alternativa no futuro das comunicacoes moveis.
O processo de construcao dos Codigos de Controle da Paridade CPq(n), usandoCodigos de Goppa, permite colocar esta importante classe de codigos algebricosclassicos no contexto da teoria moderna de codigos, pois eles tem se mostrado sereficientes.
A Teoria de Codigos Detectores e Corretores de Erros teve suas origens em 1948com os trabalhos de C. E. Shannon, R. Hamming, M. Golay e outros. No inıcioestes codigos foram criados usando unicamente conceitos da Algebra e Teoria dos
∗Email: [email protected]
09
Numeros. Posteriormente, em 1977, V. D. Goppa introduziu uma nova forma deconstruir codigos lineares usando Curvas Algebricas definidas sobre corpos finitos,hoje conhecidos como os Codigos Geometricos de Goppa.
Neste trabalho usamos algumas das ideias de Goppa para construir os Codigosde Controle da Paridade CPq(n), usando a tecnica de restricao de um codigo linear(Codigo de Goppa Racional). Em geral, se C e um codigo linear definido sobre umcorpo finito Fq, a restricao de C a Fq e dada por C|Fq = C ∩Fn
q . Apresentamos aquidois casos concretos e exibimos um exemplo em cada caso.
2 Codigos de Goppa
Faremos um breve resumo de conceitos e resultados sobre Curvas Algebricas e Cor-pos de Funcoes Algebricas. Neste contexto apresentamos um classico resultadoda Geometria Algebrica, o Teorema de Riemann-Roch, o qual permite estimar osparametros para os Codigos de Goppa.
Seja F um corpo extensao de um corpo K e se denota por F/K. A Teoria deCorpos de Funcoes Algebricas garante que se F/K e um corpo de funcoes de umavariavel, entao existe uma curva projetiva nao-singular X (unica salvo isomorfismo)cujo corpo de funcoes K(X) e isomorfo a F . Este isomorfismo permite ”traduzir”definicoes e resultados dos corpos de funcoes algebricas para curvas algebricas evice-versa.
Agora, considere uma curva algebrica nao-singular X, definida sobre o corpofinito Fq e seja F (X) seu corpo de funcoes (sobre Fq) associado. Sejam P1, · · · , Pn
lugares distintos de grau 1 em F/Fq (equivalentemente, P1, · · · , Pn pontos racionaisda curva). Considere o divisor (soma formal de lugares ou pontos da curva)
D =
n∑i=1
Pi
e seja G qualquer outro divisor de F/Fq tal que Pi /∈ supp(G), para todo i = 1, · · · , n(ou seja supp(G) ∩ supp(D) = ∅). O codigo geometrico de Goppa, associado aosdivisores D e G, esta definido por
CL(D,G) = (x(P1), · · · , x(Pn)) : x ∈ L(G) ⊆ Fnq ,
sendo L(G) um espaco vetorial (espaco de funcoes algebricas) dado por
L(G) = x ∈ F/Fq − 0 : div(x) +G ≥ 0 ∪ 0,
onde div(x) e o divisor (principal) associado a funcao (algebrica) x. O espaco L(G)e finito dimensional sobre Fq e a sua dimensao denota-se por ℓ(G). Por definicaotemos que dim(G) = dim L(G).
Observacao 2.1 Para um lugar P de grau 1 e um elemento x ∈ F , com vP (x) ≥ 0,temos que x(P ) e o valor de x em P (ou seja, x(P ) ∈ Fq e vP (x− x(P )) > 0). vPe dita a valoracao de F no lugar P .
Considerando a aplicacao
ϕ : L(G) −→ Fnq , dada por ϕ(x) = (x(P1), x(P2), · · · , x(Pn)),
observamos que a imagem de ϕ e um subespaco vetorial de Fnq . Este subespaco e
precisamente o Codigo Geometrico de Goppa. Dado que a aplicacao ϕ e Fq-linear eCL(D,G) = ϕ(L(G)), temos que CL(D,G) e um q-ario codigo linear.
10
Um importante resultado da Geometria Algebrica e o Teorema de Riemann-Roch, que apresentamos a seguir para o caso especıfico do corpo de funcoes F/Fq.
W e dito divisor canonico se deg(W ) = 2g − 2 e dim(W ) = dim L(W ) = g,onde g e o genero do corpo F/Fq (o genero g e um invariante do corpo de funcoesF/Fq ou da curva algebrica associada X).
Teorema 2.1 (Riemann-Roch) Seja W um divisor canonico de F/Fq. Entao, paraqualquer divisor A, temos que
dim(A) = deg(A) + 1− g + dim(W −A).
Pelo teorema de Riemann-Roch, para o divisorG da definicao de Codigo Geometricode Goppa temos
ℓ(G) ≥ deg(G) + 1− g,
e a igualdade vale se deg(G) ≥ 2g − 1.
O teorema de Riemann-Roch permite estimar os parametros para os CodigosGeometricos de Goppa.
Teorema 2.2 Sejam F/Fq um corpo de funcoes algebricas de genero g, e D ⊂F (Fq) (F (Fq) e o conjunto de lugares de grau 1), com |D| = n (cardinalidade deD). Seja G um divisor com g ≤ deg(G) < n e supp(G) ∩D = ∅. Entao CL(D,G)e um [n, k, d]-codigo linear sobre Fq satisfazendo:
k ≥ deg(G)− g + 1, d ≥ n− deg(G).
Mais ainda, se deg(G) ≥ 2g − 1, entao k = deg(G)− g + 1.
Se n e muito mais grande do que deg(G), entao ϕ e um mergulho em Fnq e a
dimensao k de CL(D,G) e igual a ℓ(G).
Um codigo de Goppa associado aos divisores do corpo de funcoes algebricasFq(z)/Fq e chamado de Codigo de Goppa Racional, onde z e um elemento transcen-dente de Fq. De maneira mais formal, considere um divisor D = P1+P2+ · · ·+Pn,onde os Pi sao lugares de grau 1 (pontos racionais da curva correspondente), e umoutro divisor E, de grau positivo, tal que supp(E) ∩ supp(D) = ∅. Seja f umafuncao racional no corpo das funcoes racionais Fq(z). Para cada f ∈ L(E) − 0temos que f(Pi) esta bem definido e e um elemento do corpo Fq.
Por meio do divisor D, a funcao avaliacao ϕ : L(E) −→ Fq, dada por ϕ(f) =(f(P1), · · · , f(Pn)), e linear sobre Fq e sua imagem e o subespaco linear de Fn
q , dadopor
C(D,E) = (f(P1), · · · , f(Pn)) : f ∈ L(E).
Portanto, C(D,E) e um codigo linear, chamado Codigo de Goppa Racional associ-ado aos divisores D e E.
Os Codigos de Goppa (classicos) foram introduzidos por meio de relacoes po-linomiais, generalizando os Codigos BCH (Bose-Chauhuri-Hocquenghem) (codigoscıclicos junto a um conjunto conveniente de raızes n-esimas da unidade e com umacota inferior mınima).
11
3 Codigos Controle da Paridade
Um codigo linear C e um subespaco vetorial do espaco Kn, sendo K = Fq um corpofinito de cardinalidade q e n ∈ N. Consideremos o corpo finito Fq como sendo oalfabeto do codigo C.
Todo codigo linear C possui tres parametros fundamentais, n,m e d, sendo nseu comprimento, m sua dimensao e d sua distancia mınima. Assim, dizemos queC e um [n,m, d]-codigo.
O Codigo de Controle da Paridade, CPq(n) e um subespaco vetorial de Fnq ,
obtido ao acrescentar a cada elemento do espaco Fn−1q uma ultima componente, de
modo que a soma de todas as entradas seja zero em Fq. Formalmente, CPq(n) edado por
CPq(n) = (c1, · · · , cn) : (c1, · · · , cn−1) ∈ Fn−1q ,
n∑i=1
ci = 0,
e e um [n, n− 1, 2]-codigo q-ario. Isto significa que CPq(n) e um codigo de compri-mento n, dimensao n − 1 e distancia mınima 2 e, portanto, trata-se de um codigodetector de um erro simples.
O Codigo de Controle da Paridade CPq(n) sera obtido como a restricao de umCodigo de Goppa Racional. Serao considerados dois casos, dependendo de se ocomprimento n do codigo e divisıvel ou nao pela caracterıstica do corpo finito Fq.
Teorema 3.1 Se char(Fq) - n, entao CL(D,G)|Fq = CPq(n).
Prova: Seja m um inteiro tal que n|qm−1 e seja β ∈ Fqm uma n-esima raizprimitiva da unidade. Consideremos o corpo de funcoes racionais Fqm(z)/Fqm edenotemos por Pi os zeros da funcao z − βi−1, para i = 1, · · · , n. Tambem conside-remos os divisores
G = (n− 1)P∞ − P0 e D =n∑
i=1
Pi,
onde P0 e o lugar zero e P∞ o lugar infinito de z em Fqm(z)/Fqm .
Dado que deg(G) = n− 2, o Teorema de Riemann-Roch garante que dim(G) =n − 1. Assim o conjunto z, z2, · · · , zn−1 e uma base para o espaco L(G). Agoraconsidere a aplicacao avaliacao ϕ da secao anterior. Para cada j = 1, · · · , n − 1,temos que
ϕ(zj) = (1, βj , · · · , (βn−1)j),
e como β e uma n-esima raiz primitiva da unidade, obtemos
n−1∑i=0
(βi)j = 0,
ou seja, a soma das componentes e igual a zero em Fq. Logo, todo elemento x ∈ L(G)e da forma
x =
n−1∑k=1
akzk,
12
com ak ∈ Fqm , para k = 1, · · · , n− 1. Portanto
ϕ(x) = (x(P1), x(P2), · · · , x(Pn))
=
(n−1∑k=1
akzk(P1),
n−1∑k=1
akzk(P2), · · · ,
n−1∑k=1
akzk(Pn)
)
=
(n−1∑k=1
ak,n−1∑k=1
akβk, · · · ,
n−1∑k=1
ak(βn−1)k
).
Em consequencia temos que CL(D,G) = CPqm(n), pois
n−1∑k=1
ak + · · ·+n−1∑k=1
ak(βn−1)k =
n−1∑k=1
ak
(n−1∑i=0
(βi)k)= 0,
e os vetores vj = (1, βj , · · · , (βn−2)j), para j = 1, 2, · · · , n, formam uma basepara o espaco Fn−1
qm . Assim resulta que
CL(D,G)|Fq = CPq(n).
Teorema 3.2 Se char(Fq)|n, entao CL(D,G)|Fq = CPq(n).
Prova: De maneira similar, seja m um inteiro tal que (n− 1)|qm−1 e β ∈ Fqm
uma (n−1)-esima raiz primitiva da unidade. Considere o corpo de funcoes racionaisFqm(z)/Fqm e denotemos por Pi os zeros da funcao z − βi−1, para i = 1, · · · , n− 1e por Pn := P0 o zero de z.
Consideremos os divisores
G = (n− 2)P∞ e D =n−1∑i=1
Pi,
no corpo de funcoes racionais Fqm(z)/Fqm , como no caso 1.
Como deg(G) = n−2, entao dim(G) = n−1, e portanto o conjunto 1, z, z2, · · · , zn−2e uma base para o espaco L(G). Alem disso, temos que
ϕ(1) = (1, 1, · · · , 1) e ϕ(zj) = (1, βj , · · · , (βn−2)j , 0).
Entao e claro que a soma das componentes de ϕ(1) e ϕ(zj), para j = 1, 2, · · · , n−2,e igual a zero em Fq. Desta forma, como no caso anterior, obtemos
CL(D,G)|Fq = CPq(n).
4 Aplicacoes
1) O codigo CP2(3):
Dado que 2 - 3, seja β ∈ F4 uma raiz primitiva terceira da unidade. Considereo corpo de funcoes racionais F4(z)/F4. Dado que F4 = 0, 1, β, β2, denotemos
13
os lugares de grau 1, em F4(z)/F4, por P0, P1, Pβ, Pβ2 e P∞. Agora considere osdivisores
D = P1 + Pβ + Pβ2 e G = 2P∞ − P0,
em F4(z)/F4. Assim temos que z, z2 e uma base para o espaco L(G). Segue que,dos 16 elementos de L(G), os unicos cuja imagem a respeito de ϕ estao em F3
2 sao
0, z + z2, βz + β2z2, β2z + βz2.
Desta forma, sua imagem e exatamente (0, 0, 0), (0, 1, 1), (1, 1, 0), (1, 0, 1). Assimobtemos que
CL(D,G)|F2 = CP2(3).
2) O codigo CP2(4):
Dado que 2|4, considere β ∈ F4 e o corpo de funcoes racionais F4(z)/F4, comono caso anterior. Agora considere os divisores
D = P1 + Pβ + Pβ2 + P0 e G = 2P∞,
em F4(z)/F4. Assim temos que 1, z, z2 e uma base para o espaco L(G). Segueque os unicos elementos, cuja imagem a respeito de ϕ estao em F4
2 e suas respectivasimagens, sao:
0 (0,0,0,0)
1 (1,1,1,1)
z + z2 (0,1,1,0)
βz + β2z2 (1,1,0,0)
β2z + βz2 (1,0,1,0)
z + z2 + 1 (1,0,0,1)
βz + β2z2 + 1 (0,0,1,1)
β2z + βz2 + 1 (0,1,0,1)
Obtemos assim CL(D,G)|F2 = CP2(4).
Referencias
[1] H. Stichtenoth; Algebraic Function Field and Codes, Springer-Verlag, Berlin-Heidelberg, 1993.
[2] V. D. Goppa; Geometry and Codes, Kluwer Academic Publisher, Boston 1988.
[3] J. Van Lint ; Introduction to Coding Theory, Second edition, Springer-Verlag,1992.
[4] A. Hefez e M. L. T. Vilela; Codigos Corretores de Erros, Serie de Computacaoe Matematica, IMPA, Rio de Janeiro, 2002.
[5] S. Roman; Coding and Information Theory, Springer-Verlag, N.Y., 1991.
14
Calculo Fracionario Aplicado ao Problema da
Tautocrona∗
Pedro Felipe Pavanelo Ramos e Rubens de Figueiredo Camargo.†
1 de novembro de 2012
Resumo
Apresentamos neste trabalho um estudo introdutorio sobre integrais e derivadasde ordens arbitrarias, o assim chamado calculo de ordem nao-inteira, popularizadocom o nome de Calculo Fracionario. Em particular, discutimos e resolvemos oclassico problema da tautocrona, isto e, o problema de se determinar a curva na qualo tempo gasto por um objeto para deslizar, sem friccao, em gravidade uniforme ateseu ponto de mınimo e independente de seu ponto de partida, utilizando integrais ederivadas de ordem nao inteira.
Palavras Chave: Calculo Fracionario, Integral Fracionaria, Derivada Fracionariade Riemann-Liouville.
Introducao
O calculo fracionario tem sua origem em 30 de Setembro de 1695 em uma cartaescrita por l’Hospital ao seu amigo Leibniz [13, 18], na qual o significado de umaderivada de ordem meio e proposto e discutido.
A resposta de Leibniz ao seu amigo, somada a contribuicao de inumeros bri-lhantes matematicos como Euler, Lagrange, Laplace, Fourier, Abel, Heaviside, Liou-ville, entre outros, levaram as primeiras definicoes de derivadas e integrais de ordensnao-inteiras, que no final do seculo XIX, devido primordialmente as definicoes pro-postas por Riemann-Liouville e Grunwald-Letnikov, pareciam estar completas[16].
Ate proximo ao final do seculo passado o desenvolvimento do calculo fracionariodeu-se estritamente no campo da matematica pura, sem grandes aplicacoes em ou-tras areas. Contudo, em 1969 M. Caputo, em seu livro Elasticita e Dissipazione [13],resolveu problemas de viscoelasticidade utilizando uma nova definicao, proposta porele, para a derivada de ordem fracionaria. Alem disso, o autor utilizou sua definicaopara descrever problemas de sismologia. Por outro lado, a assim chamada deri-vada fracionaria de Grunwald-Letnikov, mostrou-se bastante eficiente para resolverproblemas numericos.
A partir das definicoes introduzidas por Caputo e Grunwald-Letnikov, que temum carater local, nas ultimas decadas diversos autores mostraram que a modelagem
∗Trabalho realizado como parte do relatorio final do projeto de iniciacao cientıfica entitulado: DasTransformadas Integrais ao Calculo Fracionario aplicado a Engenharia de Producao, Bolsa Fapesp, pro-cesso 2011/07241-8. sob a Orientacao do Prof. Dr. Rubens de Figueiredo Camargo.
†Email: pedrof [email protected] e [email protected], respectivamente, curso de Engenharia deProducao, UNESP-Bauru e docente do departamento de Matematica da UNESP de Bauru.
15
feita a partir do calculo fracionario oferece uma descricao mais fina de fenomenosnaturais que aquela feita a partir do calculo usual, tendo em vista que derivadasfracionarias proporcionam uma excelente descricao para efeitos de memoria e pro-priedades hereditarias [13].
A maneira canonica de se utilizar esta importante ferramenta e substituir aderivada de ordem inteira de uma equacao diferencial parcial, que descreve umdeterminado fenomeno, por uma derivada de ordem nao-inteira. Varios resultadosimportantes e generalizacoes foram obtidos atraves desta tecnica, em diversas areasdo conhecimento, tais como: mecanica dos fluidos, fenomenos de transporte, redeseletricas, probabilidade, biomatematica, dentre outros [1, 5, 8, 9, 10, 11, 12, 13, 14,17, 18].
No presente trabalho, introduzimos os conceitos de integral e de derivada deordens nao inteira, segundo Riemann Liouville e utilizamos estes conceitos pararesolver o problema da tautocrona, que consiste em se determinar a curva na qualo tempo gasto por um objeto para deslizar, sem friccao, em gravidade uniforme ateseu ponto de mınimo e independente de seu ponto de partida [15].
1 Integral Fracionaria
Vamos introduzir o conceito de integral fracionaria como uma generalizacao daintegral de ordem inteira. Para tanto, vamos mostrar que uma integral de ordemn, de uma funcao integravel f(x), pode ser vista como um produto de convolucaode Laplace entre f(x) e a funcao Gel’fand-Shilov de ordem n. Em seguida, vamosutilizar a generalizacao do conceito de fatorial atraves da funcao gama para gene-ralizar o conceito de integral de ordem inteira.
Definicao 2: Funcao de Gel’fand-Shilov
Sejam n ∈ N e ν 6∈ N, definimos a funcao de Gel’fand-Shilov como
φn(t) :=
tn−1
(n − 1)!se t ≥ 0
0 se t < 0e φν(t) :=
tν−1
Γ(ν)se t ≥ 0
0 se t < 0.
Definicao 3: Operador Integral de ordem nDefinimos a integral de ordem inteira atraves do operador I definido como
If(t) =
∫ t
0
f(t1) dt1.
Desta forma temos
I2f(t) = I[If(t)] =
∫ t
0
∫ t1
0
f(t2) dt2 dt1.
De maneira analoga para o operador integral de ordem n, isto e,
Inf(t) =
∫ t
0
∫ t1
0
∫ t2
0
· · ·∫ tn−1
0
f(tn) dtn dtn1. . . dt2 dt1.
Definicao 4: Convlucao de Laplace
Definimos a convolucao, ou produto de convolucao de Laplace, de f(t) e g(t),denotada por f(t) ∗ g(t), como sendo [7]
16
f(t) ∗ g(t) =∫ t
0
f(t− τ)g(τ) dτ =
∫ t
0
f(τ)g(t− τ) dτ .
Teorema 1: Integral de ordem nSeja f(x) uma funcao integravel, entao
Inf(t) = φn(t) ∗ f(t) =∫ t
0
φn(t− τ)f(τ) dτ =
∫ t
0
(t− τ)n−1
(n − 1)!f(τ) dτ. (1.0.1)
Demonstracao: Vamos provar o teorema por inducao em n. Para n = 1 temosque
If(t) =
∫ t
0
f(τ) dτ =
∫ t
0
(t− τ)1−1
(1− 1)!f(τ) dτ = φ1(t) ∗ f(t).
Para concluir basta mostrar que se Inf(t) = φn(t) ∗ f(t) entao In+1f(t) =I[Inf(t)] = φn+1(t) ∗ f(t). Por hipotese indutiva temos que
In+1f(t) = I[φn(t) ∗ f(t)] =∫ t
0
φn(u) ∗ f(u)du =
∫ t
0
∫ u
0
(u− τ)n−1
(n− 1)!f(τ) dτdu.
Pelo teorema de Goursat [13] podemos mudar a ordem de integracao, i.e.,
In+1f(t) =
∫ t
0
[∫ t
τ
(u− τ)n−1
(n− 1)!f(τ) du
]
dτ.
Calculando a integral entre colchetes podemos escrever
In+1f(t) =
∫ t
0
(t− τ)n
n!f(τ) dτ = φn+1(t) ∗ f(t),
que e justamente o resultado desejado, a partir do qual apresentamos a definicaoque se segue.
Definicao 5: Integral de ordem arbitraria segundo Riemann-Liouville.
Seja f(t) uma funcao integravel, utilizamos a generalizacao da funcao fatorialpela funcao gama [4] para definirmos a integral de ordem ν de f(t), denotada porIνf(t), por
Iνf(t) = φν(t) ∗ f(t) =∫ t
0
(t− τ)ν−1
Γ(ν)f(τ) dτ. (1.0.2)
2 Derivada Fracionaria
Ha varias formas de se introduzir a derivada de ordem nao-inteira como umageneralizacao para a derivada de ordem inteira, dentre elas podemos citar a definicaode Riemann-Liouville, que e a mais conhecida e a de Caputo, que e mais restritiva,mas parece ser mais adequada para o estudo de problemas fısicos [9, 10, 11, 12]. Alemdisso, destacamos a definicao de Grunwald-Letnikov que e mais apropriada para seutilizar em problemas numericos. Estas e outras definicoes podem ser encontradasem detalhes no livro de Podlubny [18].
No presente trabalho, estamos interessados na resolucao de uma equacao inte-grodiferencial fracionaria, relacionada ao problema da tautocrona, por esta razaoapresentamos apenas a derivada fracionaria segundo Riemann-Liouville, que comoveremos, baseia-se no fato da derivacao ser o operador inverso a esquerda da inte-gracao e na lei dos expoentes de Lagrange [13].
Definicao 6: Derivada fracionaria segundo Riemann-Liouville.
Sejam f(t) uma funcao diferenciavel, m ∈ IN e α 6∈ IN tais que m− 1 < Re(α) < m.
17
A derivada de ordem α no sentido de Riemann-Liouville e definida como sendoa derivada de ordem inteira de uma integral fracionaria, de forma que a lei dosexpoentes faca sentido, isto e,
Dαf(t) = DmIm−αf(t) = Dm[φm−α(t) ∗ f(t)]. (2.0.3)
3 O Problema da Tautocrona
Apresentamos aqui a primeira aplicacao do calculo fracionario presente na li-teratura, ou seja, a solucao proposta por Abel para o problema da Tatocrona oucurva isocronica, isto e, o problema de se determinar a curva na qual o tempo gastopor um objeto para deslizar, sem friccao, em gravidade uniforme ate seu ponto demınimo e independente de seu ponto de partida.
Vamos resolver este problema de duas maneiras distintas a primeira utilizandoo calculo usual e a segunda utilizando o calculo fracionario.
A solucao proposta por Abel parte do Princıpio da Conservacao de Energiaque estabelece que a quantidade total de energia em um sistema isolado permanececonstante, ou ainda, que a soma entre energia potencial gravitacional e a energiacinetica e constante.
Por hipotese, no problema da tautocrona, a partıcula se move sem atrito e con-sequentemente sua energia cinetica e exatamente igual a diferenca entre a energiapotencial em seu ponto inicial e a energia potencial no ponto em que se encontra.Sendo m a massa do objeto, v(t) sua velocidade no instante t, y0 a altura em quefoi abandonado e y(t) a altura no instante t, sabemos que as energias cinetica epotencial sao dadas, respectivamente, por
1
2mv2 e mg y.
Como a partıcula esta restrita a mover-se sobre a curva, sua velocidade e simples-mente v = ds/dt, onde s e a distancia medida ao longo da curva. Utilizando oprincıpio da conservacao de energia podemos escrever
1
2m
(
ds
dt
)2
= mg(y0 − y)
ds
dt= ±
√
2g(y0 − y)
dt = ± ds√
2g(y0 − y)
dt = ± 1√2g
(y − y0)−1/2
(
ds
dy
)
dy,
onde a ultima passagem e devida a regra da cadeia. Note que a funcao s(y) descrevea distancia remanescente na curva em termos da altura remanecente y, como adistancia e a altura diminuem a medida que o tempo passa devemos tomar apenaso sinal negativo na equacao anterior, isto e;
dt = − 1√2g
(y − y0)−1/2
(
ds
dy
)
dy. (3.0.4)
Integrando a equacao (3.0.4) em ambos os lados de y0 a zero temos
τ = T (y0) =
∫
0
y0
dt = − 1√2g
∫
0
y0
(y − y0)−1/2
(
ds
dy
)
dy,
18
ou seja,
τ =1√2g
∫ y0
0
(y − y0)−1/2
(
ds
dy
)
dy, (3.0.5)
na qual τ e o tempo de descida. Esta equacao integrodiferencial pode ser resolvidade varias formas, a seguir apresentamos a classica solucao que utiliza a transfor-mada de Laplace e posteriormente a solucao proposta por Abel utilizando o calculofracionario.
3.1 Transformada de Laplace
Aplicando a transformada de Laplace em ambos os lados da equacao (3.0.5)podemos escrever
L[τ ] =1√2g
L
[∫ y0
0
(y0 − y)−1/2
(
ds
dy
)
dy
]
=1√2g
L
[
1√y
]
L
[
ds
dy
]
,
a ultima igualdade e devida ao teorema da convolucao.Sendo s o parametro da transformada de Laplace, temos que L[1] = 1/s e
L
[
1√y
]
=√π s−1/2, logo
L
[
ds
dy
]
=τ√2g√π
s−1/2.
A fim de explicitar a solucao, isto e, ds/dy, aplicamos a transformada de Laplaceinversa de onde obtemos
ds
dy=
τ√2g
π
1√y. (3.1.1)
3.2 Solucao de Abel
A solucao proposta por Abel baseia-se na observacao de que a integral da equacao(3.0.5), exceto pelo fator multiplicativo 1/Γ(1/2) = 1/
√π e exatamente a definicao
de integracao fracionaria de ordem 1/2. Desta forma aplicando a derivada segundoRiemann-Liouville de ordem 1/2 em ambos os lados da equacao temos1
d1/2
dy1/2τ =
√π√2g
d1/2
dy1/2
[
1
Γ(1/2)
∫ y0
0
(y − y0)−1/2
(
ds
dy
)
dy
]
. (3.2.1)
Lembrando que a derivada de Riemann-Liouvile de ordem 1/2 de 1 = y0 vale(√π y)−1/2 podemos escrever
τ√π y
=
√π√2g
(
ds
dy
)
,
ou seja,ds
dy=
τ√2g
π
1√y
1Lembre-se que a derivada fracionaria de Riemann-Liouville e o operador inverso a esquerda da integralfracionaria.
19
que e o mesmo resultado obtido na equacao (3.1.1). Note que a solucao via calculofracionario e bem mais simples e direta que a solucao classica. Essa equacao(variaveis separaveis) nos mostra que
ds =τ√2g
π
1√ydy,
de onde segue, apos a integracao, que
s(y) =2 τ
√2g
πy1/2.
3.3 Representacao grafica
A fim de ilustrar a solucao do problema da tautocrona apresentamos a seguiruma representacao grafica da curva correspondente com cinco massas abandonadasem pontos diferentes, e verificamos que todas atingem o vertice ao mesmo tempo.
0.5 1.0 1.5
-1.0
-0.8
-0.6
-0.4
-0.2
0.5 1.0 1.5
-1.0
-0.8
-0.6
-0.4
-0.2
0.5 1.0 1.5
-1.0
-0.8
-0.6
-0.4
-0.2
t = 0
t = 0.6τ
t = 0.9τ
S(t)
S(t)
S(t)
t
t
t
Figura 1: Posicao das partıculas para diferentes valores de t.
20
0.2 0.4 0.6 0.8 1.0
0.5
1.0
1.5
2.0
t
S(t)
y0 = 0.2
y0 = 0.4
y0 = 0.6
y0 = 0.8y0 = 1
Figura 2: Caminho percorrido por cada partıcula.
4 Conclusoes
Concluımos que alem de ser uma ferramente precisa para refinar a descricao defenomenos naturais, o calculo fracionario tambem pode ser utilizado para obter demaneira simples e elegante a solucao de problemas de calculo variacional, tais comoo da Tautocrona.
5 Agradecimentos
Gostarıamos de agradecer a FAPESP por ter financiado este projeto, processo2011/07241-8. Alem disso, agradecemos ao Prof. Dr. Alexys Bruno Alfonso porsua indispensavel ajuda com a parte grafica e aos pareceristas por suas valiosascontribuicoes e correcoes.
Referencias
[1] A. Luisa Soubhia, R. Figueiredo Camargo, E. Capelas de Oliveira, J. Vaz Theo-rem for Series in Thee-Parameter Mittag-Leffler Function, Fractional Calculusand Applied Analysis, bf13, (2010).
[2] H. Bernhard The Bernoulli family, In: Wussing, H. Arnold, W. Biographienbedeutender Mathematiker. Berlin, 1983.
[3] W. E. Boyce and R. C. DiPrima, Equacoes Equacoes Diferenciais Elementarese Problemas de Valores de Contorno, Oitava Edicao, LTC, Rio de Janeiro,(2006).
[4] E. Capelas de Oliveira, Funcoes Especiais com Aplicacoes, Editora Livraria daFısica, Sao Paulo, (2005).
[5] L. Debnath, Recent Applications of Fractional Calculus to Science and Engi-neering, Int. J. Math. Math. Sci., 54, 3413-3442, (2003).
[6] Felix Silva Costa ; Capelas de Oliveira, E. ; R. Figueiredo Camargo,; J.Vaz. Integrais de Mellin-Barnes e a Funcao de Fox. TEMA. Tendenciasem Matematica Aplicada e Computacional, v. 11, p. 01/01, (2011).[doi:10.5540/tema.2011.012.02.0157]
21
[7] R. Figueiredo Camargo, Do Teorema de Cauchy ao Metodo de Cagniard, Dis-sertacao de Mestrado, Universidade Estadual de Campinas, Campinas, (2005).
[8] R. Figueiredo Camargo, E. Capelas de Oliveira and F. A. M. Gomes, TheReplacement of Lotka-Volterra Model by a Formulation Involving FractionalDerivatives, R. P. 10/07, (2007).
[9] R. Figueiredo Camargo, Ary O. Chiacchio and E. Capelas de Oliveira, Differen-tiation to Fractional Orders and the Fractional Telegraph Equation, J. Math.Phys., 49, 033505, (2008) [doi:10.1063/1.2890375].
[10] R. Figueiredo Camargo, R. Charnet and E. Capelas de Oliveira, OnSome Fractional Green’s Functions, J. Math. Phys. 50, 043514 (2009).[doi:10.1063/1.3119484].
[11] R. Figueiredo Camargo, Ary O. Chiacchio, R. Charnet and E. Capelas deOliveira, Solution of the fractional Langevin equation and the Mittag-Lefflerfunctions, J. Math. Phys. 50, 063507 (2009). [doi:10.1063/1.3152608].
[12] R. Figueiredo Camargo, E. Capelas de Oliveira and J. Vaz Jr. On anoma-lous diffusion and the fractional generalized Langevin equation for a harmonicoscillator J. Math. Phys. 50, 123518 (2009). [doi:10.1063/1.3269587].
[13] R. Figueiredo Camargo, Calculo Fracionario e Aplicacoes, Tese de Doutorado,Unicamp, Campinas, (2009).
[14] R. Figueiredo Camargo,; Oliveira, Edmundo Capelas ; Vaz, Jayme. On theGeneralized Mittag-Leffler Function and its Application in a Fractional Tele-graph Equation. Mathematical Physics, Analysis and Geometry , v. 14, p. 274,(2011). [doi:10.1007/s11040-011-9100-8]
[15] R. Figueiredo Camargo e E. Capelas de Oliveira. Introducao ao Calculo deOrdem nao Inteira., livro em fase final de editoracao.
[16] K. S. Miller and B. Ross, An Introduction to the Fractional Calculus and Frac-tional Differential Equations, John Wiley & Sons, Wiley-Interscience, (1993).
[17] G. Pedroso Pires Abreu e R. Figueiredo Camargo, Modelagem Matematica: DasTransformadas Integrais ao Calculo de Ordem nao-Inteira, XXIII Congresso deIniciacao Cientıfica da UNESP, Bauru, (2011).
[18] I. Podlubny, Fractional Differential Equations, Mathematics in Science andEngineering, Vol.198, Academic Press, San Diego, (1999).
22
Metodos Iterativos para a Solucao da Equacao de Poisson
Valdirene da Rosa Rocho, Dagoberto Adriano Rizzotto Justo,Programa de Pos-Graduacao em Matematica Aplicada, PPGMap, UFRGS,
91509-900, Porto Alegre, RSE-mail: [email protected], [email protected]
Palavras-chave: Equacao de Poisson, Condicoes de Contorno de Neumann, Metodos Iterati-vos, Convergencia.
Resumo: Para aproximar a solucao da equacao de Poisson atraves do metodo de diferencasfinitas precisamos resolver um sistema linear, que pode ser resolvido atraves de um metodo ite-rativo. Para analisar a convergencia de tais metodos podemos estudar os autovalores do sistemaobtido, onde desejamos que o modulo do maior autovalor seja menor ou igual a um. Apresen-tamos neste trabalho formulas para todos os autovalores obtidos utilizando condicao de contornode Neumann. Para este problema obtem-se condicoes para que o mesmo tenha solucao baseadona integral do termo fonte.
Introducao
A equacao de Poisson e uma equacao elıptica de derivadas parciais com uma ampla utilidadeem Dinamica de Fluidos. Para resolver esta equacao pode-se utilizar varios metodos como,por exemplo, uma funcao de Green ou metodos numericos. E ainda precisa-se de adequadascondicoes de contorno.
Uma maneira de aproximar a solucao da equacao de Poisson e utilizar o metodo de diferencasfinitas, de onde obtem-se um sistema linear A~p = ~b esparso e tridiagonal no caso unidimensional.Esse sistema pode ser resolvido atraves de um metodo iterativo como Jacobi ou Gauss-Seidel,onde uma matriz de iteracao G e obtida [1, 6]. Para analisar a convergencia de tais metodospodemos estudar os autovalores da matriz G, onde desejamos que o modulo do maior autovalorseja menor ou igual a um.
O problema aliado a condicoes de contorno de Dirichlet e amplamente estudado na literaturacomo em [8, 7, 6] onde sao apresentados para alguns casos formulas para as matrizes de iteracao.Entretanto, para o problema com condicoes de Neumann nao sao encontradas formulas oucondicoes clara para a solucao de tais problemas. Apresentamos neste trabalho formulas paratodos os autovalores obtidos utilizando condicoes de contorno de Neumann.
Utilizando o produto de Kronecker pode-se relacionar as matrizes de discretizacao dos pro-blemas em uma e duas dimensoes. O espectro da matriz A e util, por exemplo, para determinara existencia de uma ou infinitas solucoes do problema discretizado. Os problemas de Neumannpossuem um autovalor igual a zero, o que nao garante solucao para o problema. Se esta existir,nao sera unica.
Ao tratar do problema de Neumann analisou-se as matrizes de iteracao a que se refereaos metodos iterativos e obteve-se que o maior autovalor e exatamente 1, e de acordo com [8]deverıamos ter que o raio espectral fosse estritamente menor do que 1 para o metodo iterativoser convergente. Seguindo [5] podemos concluir sobre quais condicoes o problema de Neumanne condicionalmente convergente, pois encontrou-se que o autovalor λ = 1 e unico, e assim seraconsiderado o segundo maior autovalor como sendo o pseudo raio espectral da matriz de iteracao.
23
1 A Equacao de Poisson
A equacao de Poisson e uma equacao de derivadas elıptica, dada por
∆p = f (x, y) ∈ Ω (1)
Segundo [2], ao utilizar-se as condicoes de contorno de Neumann, isto e,
∂p
∂n(x, y) = 0 (x, y) ∈ ∂Ω (2)
para resolver-se a equacao (1) e necessario que a hipotese de compatibilidade∫ ∫Ωf dΩ =
∫∂Ωg dS (3)
seja satisfeita para existir solucao para o problema.Ao resolver a equacao (1) para o caso unidimensional com condicoes de contorno de Neumann,
discretizamos-a porp1 − p0
h= 0, o que implica que p0 = p1, (4)
epn+1 − pn
h= 0, o que implica que pn+1 = pn. (5)
Portanto discretizando-se o problema contınuo (1) no caso unidimensional tem-se
pi+1 − 2pi + pi−1
h2= fi (6)
A partir de (6) obtem-se um sistema de equacoes lineares. Considerando as condicoes decontorno (4) e (5) este sistema admite a forma matricial
−2 2 . . . 01 −2 1
. . . . . . . . .1 −2 1
0 . . . 2 −2
.
p0
p1...pnpn+1
=
0
h2f1...
h2fn0
. (7)
Dessa forma, o problema de encontrar uma aproximacao para a funcao p(x) que satisfaz (1)se reduz em encontrar um vetor ~p = (p1, ..., pn) que satisfaz a equacao matricial (7). Em geral,sabe-se que a solucao do problema nao e unica, isso implica que se o sistema A~p = ~b tem solucao~p, entao ~p+ α.~1 tambem sera solucao desse sistema, ou seja,
A(~p+ α.~1) = A~p+ αA~1 = ~b+ 0 = ~b (8)
Para o caso bidimensional a discretizacao da equacao de Poisson atraves do metodo diferencasfinitas centradas discretizada no ponto xij , yij torna-se a formula de cinco pontos, que e dadapor
pi+1,j + pi−1,j − 4pi, j + pi,j+1 + pi,j−1
h2= fi,j . (9)
Assim, obtem-se o sistema linear (L~p = ~b). Utilizando o produto de kronecker reescreve-sea matriz do caso bidimensional [4]
L = (A⊗ I) + (I ⊗A) (10)
24
onde A e uma matriz tridiagonal, ou seja, a matriz do sistema (7). Deste modo, descreve-se
L = A⊕A =
A− 2I 2I . . . 0I A− 2I I
. . . . . . . . .I A− 2I I
0 . . . 2I A− 2I
(11)
onde o bloco I e a matriz identidade e o bloco A− 2I de (11) e dado por
A− 2I =
−4 2 . . . 01 −4 1
. . . . . . . . .1 −4 1
0 . . . 2 −4
(12)
2 Sistema Linear e Metodos Iterativos
Os metodos iterativos sao utilizados para encontrar a solucao de um sistema de equacoes linearesdescritas como
A~x = ~b (13)
com n equacoes e n incognitas. Uma simples aproximacao para uma solucao iterativa de umsistema linear e reescrever (13) como uma iteracao de ponto fixo linear [8] de modo a obter umaequacao iterativa como
~xk+1 = G~xk + ~f (14)
que convirja para a solucao do mesmo, onde G ∈ Rnxn e a matriz de iteracao do metodo iterativo(14) e ~f ∈ Rn.
O metodo de Jacobi e de Gauss-Seidel por ser splitting nos permite decompor A, assimdeterminou-se que GJ = D−1(E + F ) e GGS = (D−E)−1F , respectivamente. Seja A dada por(7) e portanto como GJ = D−1(E + F ) obtem-se que
D =
−2 . . . 0
−2. . .
−20 . . . −2
(15)
e
− (E + F ) =
0 2 . . . 01 0 1
. . . . . . . . .1 0 1
0 . . . 2 0
. (16)
Portanto GJ e
GJ =12
0 2 . . . 01 0 1
. . . . . . . . .1 0 1
0 . . . 2 0
. (17)
25
Pode-se verificar numericamente que os autovalores desta matriz de acordo com [3], sao
λk = cos[
(k − 1)πn− 1
]para 1 ≤ k ≤ n. (18)
Segundo [8], os autovalores da matriz de iteracao do metodo Gauss-Seidel sao dados por
µk =[cos(k − 1n− 1
)π
]2
para 1 ≤ k ≤⌊n
2
⌋(19)
e os demais autovalores sao iguais a zero.O teorema 1 revela que os autovalores da matriz de discretizacao bidimensional sao justa-
mente dados pela soma dos autovalores da matriz de discretizacao unidimensional.
Teorema 1 Se A ∈ Rn×n tem autovalores λi e B ∈ Rm×m tem autovalores µj, entao a somade kronecker A⊕B = (Im ⊗A) + (B ⊗ In) tem mn autovalores
λ1 + µ1, ..., λ1 + µm, λ2 + µ1, ..., λ2 + µm, ..., λn + µm (20)
Para determinar os autovalores de G no caso bidimensional usou-se o teorema 1, ja quesabe-se quem sao os autovalores de GJ no caso unidimensional. Para determina-los e necessarioconhecer quem sao os autovalores no caso unidimensional e estes sao dados por (18). Assim osautovalores sao fornecidos por
λi,j =cos (i−1)π
n−1 + cos (i−1)πn−1
2para 1 ≤ i, j ≤ n (21)
e os autovalores da G do metodo de Gauss-Seidel sao
µi,j =
[cos (i−1)π
n−1 + cos (i−1)πn−1
2
]2
para 1 ≤ i, j ≤⌊n
2
⌋(22)
e os demais sao iguais a zero.Quando k = 1 em (21) obtem-se o autovalor λ1 = 1 e quando k = n o autovalor e λn = −1
(todos os outros autovalores satisfazem |λ| < 1). Sendo este em modulo o raio espectral damatriz de iteracao para o metodo iterativo de Jacobi e tambem o de Gauss-Seidel. E assim osmetodos iterativos em questao, sao condicionalmente convergente, e o seu raio espectral e dadopelo maior autovalor diferente de |λ| = 1 [5, 9].
3 Conclusao
A matriz gerada a partir do problema de Neumann por ser uma matriz esparsa e tridiagonalpermite com facilidade encontrar os autovalores de cada uma das matrizes, a qual seguiu-se aliteratura de [9]. Apos a obtencao do sistema A~p = ~b investigou-se os metodos iterativos Jacobie Gauss-Seidel para a resolucao de sistemas lineares gerado a partir da equacao de Poisson. Estaequacao exige a utilizacao de condicoes de contorno adequadas, a escolha destas nao e unica ea convergencia do metodo depende delas. Ao tratar do problema de Neumann analisou-se asmatrizes de iteracao a que se refere aos metodos iterativos e obteve-se que o maior autovalor eexatamente 1, e de acordo com [8] deverıamos ter que o raio espectral fosse estritamente menordo que 1 para o metodo iterativo ser convergente. A literatura de [5] apresenta teoremas edefinicoes na qual conseguiu-se encontrar uma formula geral para os autovalores e entao concluı-se que os metodos usados para encontrar a solucao do problema de Neumann e condicionalmenteconvergente, pois encontrou-se que o autovalor λ = 1 e unico, e assim sera considerado o segundomaior autovalor como sendo o pseudo raio espectral da matriz de iteracao.
26
Referencias
[1] C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, (1995).
[2] H. J. van Linde, High-order Finite-Difference Methods for Poisson’s Equation, JournalMathematics of Computation, 126, (1974), 369–391.
[3] M. Neumann and R. J. Plemmons, Convergent Nonnegative Matrices and Iterative Methodsfor Consisten Linear Systems, Journal Numerical Mathematics, 31, (1978) 265–279.
[4] R. J. Plemmons, Regular Splittings and the discrete Poisson-Neumann problem, Numer.Math 25, 25 (1976) 153-161.
[5] C. Pozrikidis, A note on the regularization of the discrete Poisson-Neumann problem, Jour-nal of Computational Physics, 172 (2001) 917–923.
[6] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, 2nd ed., (2000).
[7] J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, SIAM, 2nded., (2004).
[8] D. M. Young, Iterative Solution of Large Linear Systems, Academic Press, Inc, (1988).
[9] W-C. Yueh, Eigenvalues of Several Tridiagonal Matrices, Applied Mathematics E-Notes, 5(2005) 66-74.
27
Modelo dinamico incluindo Manejo Integrado dePragas (MIP) e estrutura espacial no combate a
Diaphorina Citri∗
Priscila Azevedo da Silveira †
26 de outubro de 2012
Resumo: Propomos um modelo dinamico para o combate a Diaphorina Citri, praga quecostuma atacar culturas cıtricas em diversas regioes do Brasil, principalmente no Estadode Sao Paulo. A estrategia adotada e o Manejo Integrado de Pragas (MIP), atravesda aplicacao de pesticida e da insercao de parasitoides/predadores desta praga, levandoem consideracao a dispersao dos indivıduos adultos. Estudos analıticos e simulacoesnumericas sao apresentados, evidenciando os efeitos do MIP sobre o agroecossistema,com e sem estrutura espacial.
Palavras Chave: Manejo Integrado de Pragas, Diaphorina Citri, Modelo Matematico,Estabilidade Global.
Abstract: We propose a dynamic model for combating Diaphorina citri, a pest thatoften attack citrus crops in several regions of Brazil, especially in Sao Paulo. The strategyadopted is the Integrated Pest Management (IPM) through the application of pesticidesand release of parasitoids/predators of this pest, taking into account the dispersion ofadults. Analytical studies and numerical simulations are presented, showing the effectsof IPM on the agroecosystem, with and without spatial structure.
Keywords: Integrated Pest Management, Diaphorina Citri, Mathematical Model, GlobalStability.
Introducao
O uso de ferramentas matematicas na busca de estrategias que garantam maior efetividadeno controle de pragas em agroecossistemas tem sido objeto de inumeros estudos nasultimas decadas.
Para contornar os efeitos nocivos do uso excessivo de pesticidas agrıcolas no ambientee na saude humana, o controle biologico vem servindo como uma alternativa ao uso deprodutos quımicos no combate as pragas [2]. O MIP e uma estrategia de controle de longoprazo, que associa taticas da biologia, da quımica, bem como de cultivo, e cujo objetivonao e o de erradicar completamente a populacao de pragas, mas sim reduzi-la para nıveistoleraveis, abaixo de um Limiar Economico (LE) [10] [4].
∗Trabalho realizado como parte da tese de doutorado da autora sob Orientacao da Profa Dra. MariaCristina Varriale e Coorientacao da Profa Dra. Claudia Pio Ferreira
†Email: [email protected] - Programa de Pos-Graduacao em Matematica Aplicada - UFRGS
28
Experimentalmente a estrategia de controle atraves do MIP, alem de resultar nareducao da quantidade de inseticida a ser utilizado, tem se mostrado mais efetiva doque os metodos classicos, tais como apenas controle biologico (liberacao de predadores dapraga) ou apenas controle quımico (uso de pesticidas).
A inclusao de uma estrutura espacial no modelo matematico e, sem duvida, de fun-damental importancia, para o sucesso das estrategias de MIP, porque dependendo dastaxas e dos tipos de movimentacao (difusao, taxia, conveccao, ...) das pragas e dos seuspredadores/parasitoides, a dispersao das populacoes pode ser favoravel, ou, ao contrario,prejudicial ao controle de pragas. Alem disso, a heterogeneidade espacial pode estar pre-sente na distribuicao da populacao de pragas, sugerindo que o seu controle seja gerenciadode forma distinta, dependendo da localizacao espacial [3].
O Brasil e o maior produtor de laranjas do mundo, sendo o Estado de Sao Pauloresponsavel por 80% desta producao. Atualmente a principal praga de citros e a Di-aphorina citri (Hemiptera: Psyllidae), que transmite a bacteria Candidatus Liberibacteramericanus, causadora do “greening” - doenca atualmente presente em 18, 57% dos talhoesdo parque citrıcola de SP [5].
O parasitoide Tamarixia radiata mostra-se um bom candidato a agente de controlebiologico para a Diaphorina citri, seus ovos sao colocados dentro da ninfa do psilıdeo(em geral, um ovo por psilıdeo) e o novo parasitoide que emerge da ninfa parasitada sealimenta da mesma levando-a a morte. Alem de atuarem como parasitoides, as femeassao predadoras [6], podendo alimentar-se de ninfas da praga em fases mais jovens. Assim,somando o parasitismo a predacao, cada femea e capaz de destruir, ao longo de sua vida,uma quantidade significativa de ninfas de Diaphorina citri.
1 Formulacao e Analise do Modelo que Descreve a
Dinamica Local da D. citri e do parasitoide T.
radiata
Neste trabalho, propomos medidas para o manejo da praga Diaphorina citri em umcultivo de laranjas. Para isso, consideramos tres estagios de desenvolvimento para apraga, cujas densidades populacionais sao representadas respectivamente por: N1(t), dafase mais jovem, N2(t), da fase intermediaria, e F (t), da fase adulta. Seu inimigo natural,o inseto especialista Tamarixia radiata para o qual consideramos apenas a fase adultacom densidade populacional P (t), preda e/ou parasita, respectivamente, N1(t) e N2(t).Supondo:
1. η(1 − N1/K), a taxa de oviposicao (ovos viaveis) por femea F , onde η e a taxade oviposicao intrınseca, e K a capacidade de suporte do meio para a fase N1 -relacionada com o numero de nutrientes, espaco, etc.;
2. σN1 e σN2 , as taxas per capita com a qual N1 passa para a fase N2, e com a qual N2
passa para a fase adulta, respectivamente;
3. ψ, a proporcao dos insetos vetores adultos, que sao femeas;
4. γN1/(1+ϕN1), a taxa por predador com a qual uma presa N1 e predada, resultantede uma probabilidade N1/(1 + ϕN1) de encontro de um inseto da especie P , comum inseto vetor na fase N1;
29
5. αN2/(1 + βN2), a taxa por parasitoide femea com a qual o hospedeiro N2 e par-asitado, resultante de uma probabilidade N2/(1 + βN2) de encontro de um insetofemea da especie P , com um inseto vetor na fase N2;
6. ω e θ, os fatores de conversao de presas para predador, e de hospedeiros para para-sitoide, respectivamente;
7. µN1, µN2
, µF e µP , as taxas de mortalidade per capita de cada uma das populacoes,
escrevemos o sistema que descreve a dinamica populacional natural, isto e, semintervencao humana, sob a forma:
dN1
dt= ηF (1− N1
K)− (σN1 + µN1
)N1 − γN1P1+ϕN1
,
dN2
dt= σN1N1 − (σN2 + µN2
)N2 − αN2P1+βN2
,
dFdt
= ψσN2N2 − µFF,
dPdt
= ωγN1P1+ϕN1
+ θαN2P1+βN2
− µPP.
(1)
Neste modelo nao estamos considerando especificamente a evolucao da cultura delaranjas.
A partir deste sistema, definindo as variaveis adimensionais
τ = ηt, n1 = ϕN1, n2 = βN2, f =1
KF e p =
γ
αKP,
obtemos um sistema adimensional correspondente com os seguintes parametros adimen-sionais:
a = Kϕ, b1 =σN1
η, c1 =
µN1
η, α1 =
αKη, d = β
η,
b2 =σN2
η, c2 =
µN2
η, α2 =
αγ, α3 =
ψKβ,
ξ1 =µFη, α4 =
γωηϕ, α5 =
αθηβ
e ξ2 =µPη.
Este sistema adimensional apresenta tres pontos de equilıbrio:
E0 = (0, 0, 0, 0), E1 = (n11, n
12, f
1, 0) e E2 = (n∗1, n
∗2, f
∗, p∗),
onde n11, n
12, f
1, n∗1, n
∗2, f
∗ e p∗ sao combinacoes especıficas dos parametros do sistema.O equilıbrio E0, que corresponde a extincao de todas as especies, e biologicamente
viavel para qualquer escolha dos parametros positivos do sistema. Este equilıbrio e local-mente estavel se, e somente se,
R0 < 1,
onde
R0 =adb1b2α3
ξ1(b1 + c1)(b2 + c2). (2)
O equilıbrio E1, que corresponde a extincao da populacao de predadores e persistenciadas populacoes de pragas, e biologicamente viavel se, e somente se,
R0 > 1.
Este equilıbrio e localmente estavel se, e somente se,
R0 > 1
30
eR1 < 1,
onde R0 e dado em (2) e
R1 =Aα4[α3b2(b2 + c2) + A] + Aα5[db1b2α3 + A]
ξ2(db1b2α3 + A)[α3b2(b2 + c2) + A], (3)
comA = adb1b2α3 − ξ1(b1 + c1)(b2 + c2). (4)
As condicoes
b2α2α3α4[χ(1 + a)− aα4] < ξ1χ(b2 + c2)(α4 − χ), (5)χ+α5−α4
α4−χ < n∗2 <
ξ2−χ , se χ < 0
n∗2 >
χ+α5−α4
α4−χ , se 0 < χ < α4
, onde χ ≡ ξ2 − α5, (6)
e
n∗1 < a
α3b2n∗2
α3b2n∗2 + ξ1(b1 + c1)
, (7)
garantem que o equilıbrio E2, de coexistencia de todas as especies, seja viavel. Esteequilıbrio e localmente estavel se
ξ1(a11a22 − a14a41 − a24a42) + a11a24a42+
+a22a14a41 − db1a14a42 > (a− n∗1)db1b2α3,
(8)
a11(a11a22 − a14a41) + a22(a11a22 − a24a42)− db1a14a42 < 0 (9)
e(a− n∗
1)α3b2 < x2, (10)
onde a11, a14, a22, a24, a41, a42 e x2 sao combinacoes especıficas dos parametros adimen-sionais.
Investigamos a estabilidade global do ponto de equilıbrio de coexistencia das especies,visto que, em geral, em sistemas biologicos estamos interessados em manter o equilıbrioecologico do meio que esta sendo analisado. Atraves de uma tecnica que envolve funcoesde Lyapunov [9] [7] determinamos 2 subconjuntos da bacia de atracao do equilıbrio E2:
W1 = (n1, n2, f, p) ∈ R4+ : n1 > a, f + n1 < f ∗ + n∗
1, n2 > n∗2,
p > p∗, p− n1 < p∗ − n∗1 ∧ p− n2 < p∗ − n∗
2
eW2 = (n1, n2, f, p) ∈ R4
+ : n1 < n∗1, f > f ∗, f + n1 > f ∗ + n∗
1, n2 < n∗2,
p < p∗, p− n1 > p∗ − n∗1 ∧ p− n2 > p∗ − n∗
2.Podemos ilustrar estes resultados atraves de simulacoes numericas.Todas as simulacoes numericas deste trabalho foram obtidas a partir da implementacao
do metodo de Runge-Kutta de quarta ordem para a aproximacao das solucoes do sistema,utilizando o software Matlab.
Em todas as figuras apresentadas neste trabalho consideramos os seguintes valorespara os parametros: a = 1, α1 = 0.3, α2 = 0.5, α3 = 2 , α4 = 0.6, α5 = 0.3, b1 = 0.3,b2 = 0.1, c1 = 0.01, c2 = 0.01, d = 0.6, ξ1 = 0.3 e ξ2 = 0.4.
31
Na Figura 1, escolhendo os parametros do modelo de modo que as condicoes de ex-istencia e estabilidade local do equilıbrio E2 sejam satisfeitas e tomando a condicao inicial(n0
1, n02, f
0, p0) = (1.1, 2, 0.1, 0.3) em W1, verificamos que, de fato, as quatro populacoescoexistem evoluindo para um equilıbrio assintotico constante.
0 100 200 300 4000
0.5
1
1.5
2
Tempo
De
nsid
ad
es d
e p
rag
as e
pre
da
do
res
Figura 1: Densidades populacionais considerando a dinamica natural dos indivıduos sem estru-
tura espacial; n1 (linha pontilhada); n2 (linha tracejada); f (linha mista); p (linha contınua).
2 A Inclusao do Manejo Intregrado de Pragas (MIP)
no Modelo que Descreve a Dinamica Local da D.
citri e do parasitoide T. radiata
O MIP, que estamos propondo, consiste na aplicacao de inseticida e na liberacao deuma certa quantidade do inimigo natural da praga, quando constatado que a quantidadecrescente da densidade de femeas adultas ultrapassa LE. Desta forma, o sistema, levandoem consideracao a intervencao humana, passa a ser descrito da seguinte forma:
• Enquanto a populacao F estiver abaixo de LE, o sistema de equacoes e o sistemaque descreve a dinamica populacional natural;
• No instante τ em que a populacao de pragas atingir LE, a intervencao humanaimplicara em reduzir cada uma das populacoes atraves de uma taxa per capitaδ (dependente do tamanho da populacao de femeas adultas), e aumentar em ρ apopulacao de inimigos naturais, isto e, trabalha-se com o problema de valor inicialcomposto pelo sistema de equacoes diferenciais (1), para t ≥ τ , juntamente com asseguintes condicoes iniciais:
32
N1(τ+) = (1− δ)N1(τ),
N2(τ+) = (1− δ)N2(τ), (11)
F (τ+) = (1− δ)F (τ),
P (τ+) = (1− δ)P (τ) + ρ.
A intervencao acima sera repetida cada vez que a populacao de femeas adultas atingirLE.
A Figura 2 ilustra a aplicacao do MIP. Escolhendo, para os parametros e condicaoinicial, os mesmos valores utilizados na Figura 1 e estabelecendo LE = 30, que correspondea um valor de f = 0.3 < f ∗, verificamos que a tecnica e eficiente, mantendo a densidadede femeas adultas abaixo de LE. As densidades de ninfas de pragas no primeiro e nosegundo estagio sao mantidas abaixo do correspondente valor de equilıbrio.
0 200 400 600 8000
0.5
1
1.5
2
Tempo
De
nsid
ad
es d
e p
rag
as e
pre
da
do
res
Figura 2: Densidades populacionais considerando MIP incorporado a dinamica natural dos
indivıduos sem estrutura espacial; n1 (linha pontilhada); n2 (linha tracejada); f (linha mista);
p (linha contınua).
3 A Inclusao da Estrutura Espacial no Modelo que
Descreve a Dinamica da D. citri e do parasitoide
T. radiata
A seguir incluimos a estrutura espacial, estudando o efeito das regras de movimentacaoenvolvidas no agrossistema em questao. O habitat e representado por uma malha (ma-triz) bidimensional de tamanho L × L, com L2 sıtios (ou “patches”) discretos, cada umrepresentando uma laranjeira. Dentro de cada sıtio de coordenadas (i, j), utilizamos paraa dinamica vital o modelo que descreve a dinamica populacional natural. Assumimosque as populacoes estejam homogeneamente distribuıdas dentro de cada sıtio. A estadinamica acrescentamos regras para a movimentacao das populacoes de adultos entre os
33
diversos sıtios organizados em LM
talhoes de tamanho M ×M sıtios. As fronteiras saosupostas reflexivas. Consideramos dois tipos de movimentacao: difusao e taxia local.
Na difusao, uma fracao constante de femeas adultas (predadores), 0 < DF < 1 (0 <DP < 1), se desloca de cada sıtio e se distribui igualmente entre os oito vizinhos maisproximos (vizinhanca de Moore). A cada passo de tempo ocorrem duas atualizacoes nasdensidades dos indivıduos adultos das populacoes: uma atualizacao intermediaria daspopulacoes F e P , referente a movimentacao, que e dada pelas expressoes
Fint(i, j) = (1−DF )F (i, j) +DF
8
∑(m,n)∈V (i,j)
F (m,n), (12)
Pint(i, j) = (1−DP )P (i, j) +DP
8
∑(m,n)∈V (i,j)
P (m,n), (13)
onde V (i, j) = (i− 1, j + 1), (i, j + 1), (i+ 1, j + 1), (i+ 1, j), (i+ 1, j − 1), (i, j − 1), (i−1, j − 1), (i− 1, j).
No caso em que femeas adultas e predadores se movimentam por taxia local, umaespecie consegue perceber a densidade da outra no sıtio em que se encontra. Evidente-mente, o coeficiente de dispersao das pragas e funcao crescente da densidade de predadorese o coeficiente de dispersao dos predadores e funcao decrescente da densidade de pragas[8]. Desta forma, consideramos as funcoes
DF (P (i, j)) =P (i, j)
1 + P (i, j)(14)
e
DP (N1(i, j), N2(i, j)) =1
1 +N1(i, j) +N2(i, j)(15)
para os coeficientes de dispersao das femeas adultas da praga e dos predadores, respecti-vamente, em cada sıtio (i, j) da malha.
Quando acrescentamos ao sistema natural (sem intervencao humana) a dispersao pordifusao das femeas adultas e dos predadores construimos a Figura 3, onde apresentamos osequilıbrios que representam a soma das densidades de cada uma das populacoes em todosos sıtios do domınio. Alem disso, podemos observar que todas as populacoes tendem a sedistribuir homogeneamente pelo ambiente (caracterıstica tıpica da difusao). Neste caso,utilizamos para os parametros os valores explicitados anteriormente e para a condicaoinicial: n0
1 = 1.1, n02 = 2, f0 = 0.1 e p0 = 0.3 indivıduos em um sıtio no centro do
domınio, e todos os outros sıtios vazios.Quando aplicamos o MIP no modelo com estrutura espacial e dispersao por difusao,
podemos ver na Figura 4, que as densidades totais das populacoes de pragas passama apresentar oscilacoes com pequenas amplitudes em nıveis populacionais mais baixosdo que as densidades de equilıbrio da Figura 3. Ja a densidade total de predadoresoscila de maneira caotica. Neste caso todas as populacoes apresentam uma distribuicaoespacial heterogenea, ilustrando a interferencia do MIP na tendencia homogeneizadora dadifusao. A Figura 5 mostra a distribuicao espacial das femeas adultas apos 2000 iteracoestemporais. Para as Figuras 4 e 5 consideramos os mesmos valores dos parametros e dascondicoes iniciais da Figura 3.
No caso em que consideramos a movimentacao por taxia local o MIP tambem mostrou-se eficiente.
34
0 500 1000 1500 2000 2500 30000
1000
2000
3000
4000
Tempo
De
ns.
tota
is d
e p
rag
as e
pre
d.
Figura 3: Densidades populacionais considerando a dinamica natural dos indivıduos com dis-
persao por difusao para femeas adultas e predadores; n1 (linha pontilhada); n2 (linha tracejada);
f (linha mista); p (linha contınua).
0 200 400 600 800 1000 12000
500
1000
1500
2000
Tempo
De
ns.
to
tais
de
pra
ga
s e
pre
d.
Figura 4: Densidades populacionais aplicando MIP a dinamica natural dos indivıduos com dis-
persao por difusao para femeas adultas e predadores; n1 (linha pontilhada); n2 (linha tracejada);
f (linha mista); p (linha contınua).
Conclusao
De modo geral podemos concluir que e de fundamental importancia considerar a dis-tribuicao espacial dos indivıduos na modelagem, ja que em agroecossistemas naturais amaioria dos indivıduos adultos (pragas e/ou predadores) tem a capacidade de se movimen-tar pelo ambiente. Alem disso, os diferentes tipos de movimentacao considerados podemprovocar mudancas significativas no comportamento dos indivıduos. No modelo propostopara o caso difusivo, o MIP mostrou-se eficiente no combate a Diaphorina Citri, pois nassimulacoes podemos ver uma reducao nas densidades de pragas a nıveis toleraveis e nao
35
Solução n1 (t=2000)
10 20 30 40 50
10
20
30
40
50
Solução n2 (t=2000)
10 20 30 40 50
10
20
30
40
50
Solução f (t=2000)
10 20 30 40 50
10
20
30
40
50
Solução p (t=2000)
10 20 30 40 50
10
20
30
40
50
−0.65
−0.6
−0.55
−0.5
−0.45
−0.46
−0.455
−0.45
−0.445
−0.425
−0.42
−0.415
−0.41
−0.3
−0.299
−0.298
−0.297
Figura 5: Distribuicao espacial das populacoes apos 2000 iteracoes temporais.
a extincao destas, mantendo assim o equilıbrio ecologico. O caso da dispersao por taxialocal ainda esta sendo analisado.
Pretendemos, futuramente, aprimorar este modelo apresentando um tipo de movi-mentacao mais elaborada para os indivıduos adultos, por exemplo, considerando o fatode que estes conseguem perceber a densidade da outra especie em uma determinada viz-inhanca do sıtio onde estao. Alem disso, estamos comecando a trabalhar com medidas decontrole otimo, buscando, atraves de ferramentas matematicas, obter planejamentos maiseficientes e economicos para a aplicacao das tecnicas de controle da praga.
Referencias
[1]
[2] S. Bhattacharyya and D. K. Bhattacharya, An Improved integrated pest managementmodel under 2-control parameters (sterile male and pesticide), Math. Biosc., 209(2007) 256-281.
[3] E. A. B. F. Lima, C. P. Ferreira and W. A. C. Godoy, Ecological modeling and pestpopulation management: a possible and necessary connection in a changing world,Neotrop. Entomol. [online], 38(6) (2009) 699-707.
[4] R. F. Norris, E. P. Caswell-Chen and M. Kogan, ”Concepts in integrated pest man-agement”, Prentice Hall, New Jersey, 2003.
[5] P. E. B. Paiva, ”Distribuicao espacial e temporal, inimigos naturais e tabela de vidaecologica de Diaphorina citri Kuwayama (Hemiptera: Psyllidae) em citros em SaoPaulo”, Tese de Doutorado, Escola superior de agricultura Luiz de Queiroz-ESALQ,2009.
[6] J. A. Qureshi, M. E. Rogers, D. G. Hall and P. A. Stansly, Incidence of invasiveDiaphorina citri (Hemiptera: Psyllidae) and its introduced parasitoid tamarixia ra-
36
diata (Hymenoptera: Eulophidae) in Florida citrus, Jour. of Econ. Entomol., 215(1)(2009) 247-256.
[7] C. Robinson, ”Dynamical systems - stability, symbolic dynamics and chaos”, CRCPress, USA, 1995.
[8] P. A. Silveira, ”Perseguicao e Fuga em Modelos Presa-Predador”, Dissertacao deMestrado, UFSM, 2010.
[9] J. Sotomayor, ”Licoes de equacoes diferenciais ordinarias”, IMPA - Projeto Euclides,Rio de Janeiro, 1979.
[10] S. Tang and R. A. Cheke, Models for integrated pest control and their biologicalimplications, Math. Biosc., 215 (2008) 115-125.
37
