Rosângela dos Santos Oliveira
Aplicação do Metódo de Diferenças Finitas Em
Um Domínio Retangular
Sinop-MT, Brasil
2014
Rosângela dos Santos Oliveira
Aplicação do Metódo de Diferenças Finitas Em Um
Domínio Retangular
Trabalho de Conclusão de Curso apresentadoà Banca Examinadora do Departamento deMatemática - UNEMAT, Campus Universi-tário de Sinop, como requisito parcial paraa obtenção do título de Licenciado em Mate-mática.
Universidade do Estado de Mato Grosso – Unemat
Faculdade de Ciências Exatas
Programa de Graduação
Orientador: Miguel Tadayuki Koga
Sinop-MT, Brasil
2014
Rosângela dos Santos Oliveira
Aplicação do Metódo de Diferenças Finitas Em UmDomínio Retangular
Trabalho de Conclusão de Curso apresentadoà Banca Examinadora do Departamento deMatemática - UNEMAT, Campus Universi-tário de Sinop, como requisito parcial paraa obtenção do título de Licenciado em Mate-mática.
Trabalho aprovado. Sinop-MT, Brasil, 19 de novembro de 2014.
Miguel Tadayuki KogaOrientador
UNEMAT - Campus Universitário de Sinop
Raul Abreu de AssisProfessor Avaliador
UNEMAT - Campus Universitário de Sinop
Rodrigo Bruno ZaninProfessor Avaliador
UNEMAT - Campus Universitário de Sinop
Odacir Elias Vieira MarquesProfessor Presidente da Banca
UNEMAT - Campus Universitário de Sinop
A "toda"minha família, ao meu Amor Rodrigo incluindo sua família e aos amigos.
Agradecimentos
Agradeço a todos os professores do curso de Licenciatura em Matemática, que me
auxiliarão aos momentos de dúvidas, proporcionado experiências que contribuirão para a
minha formação intelectual, acadêmica, profissional e pessoal.
Agradeço aos professores que tive aula Celma, Dacir, Daniel, Denizalde, Falchetti,
Graci, Jean, Luciene, Marion, Nadisson, Odair que me auxiliou em algumas dúvidas sobre
o latex, Milton, Rosa Caroline, Thiélide, Vera e também aqueles que não citei os nomes.
Agradeco a Professora Luciana, por todo apoio e auxilio nas disciplinas e também sobre
uma frase que ela disse, que eu sempre lembrava em todos os momento difíceis, que no
fim das situações com esforço, dedicação e disciplina sempre da certo e ao Professor Raul
por auxiliar para os esclarecimentos das dúvidas em relação aos cálculos do método de
Crank-Nicolson.
Agradeço ao meu orientador Miguel Tadayuki Koga, por todo o auxilio, compreen-
são e esclarecimentos para construção do nosso trabalho, pois graças a ele pude conhecer
um pouco de Biomatemática. E por mostrar que é necessário ter vontade de aprender
sempre e correr atrás dos objetivos que sonhamos.
Agradeço aos professores que trabalhei no projeto Olimpíada de Matemática
da UNEMAT-Campus de Sinop, sendo eles, Chiara, Emivam, João, Galois, Miguel meu
orientador e Polyanna, por que eles me mostraram através da convivência do dia a dia,
um pouco de como é a vida de um professor os lados positivos e negativos, e também o
quanto é importante ser dedicado, ter disciplina, ter foco, força de vontade para conseguir
realizar as metas que almejamos de modo geral em nossa vida.
Agradeço a todos os amigos de longa datas que contribuíram direto e indiretamente
para a minha formação, e os colegas que convivi esses anos no curso de graduação, pelas
horas de estudos, pela ajuda mutua nos momentos de dificuldade das matérias, pelo apoio e
incentivo. Sendo eles, Alini, Amilto, Ana Caroline, Andre, Daisa e filhos, Andrei, Andressa,
Anna, Antônio Cesar, Cicero Junior, Elisângela, Fernanda, Gisele, Harry, Janaina, Jeniffer,
Joice, Karime, Karine, Kelly, Leonice, Luciana, Mayele, Maiara, Marlini, Marissa, Maxsuel,
Rafaella, Rosanise, Silvani, Tatiele Caroline, Thiago Walter e Vanessa, pois desde que
alguns deles começaram o curso os acompanhei e alguns tomaram se meus amigos.
Agradeço a toda minha "família", meu pai Clovis, minha mãe Ivone e meu Irmão
Michael e aos outros irmãos não citarei nome para que nenhum fique chateado, incluindo a
"família"do Rodrigo que é o amor da minha vida, ele uma é das pessoas mais importantes
da minha vida, e um dos motivos, e uma das razões por eu querer viver, me dedicar e
crescer a todo novo dia de forma intelectual, profissional e pessoal, enfim de maneira geral,
portanto agradeço a todos eles da minha família que amo muito por todo apoio, incentivo
auxilio em todos esses anos, por que graças a eles que me motivam e incentivavam a me
dedicar sempre aos estudos e por eles me aguentar em todas as provas e os fins de semestres
da faculdade.
"Saber muito não lhe torna inteligente. A inteligência se traduz na forma que você recolhe,
julga, maneja e, sobretudo, onde e como aplica esta informação."
(Carl Sagan)
ResumoEste trabalho visa apresentar alguns estudos desenvolvidos na área da Biomatemática
que vem sendo discutido pelo grupo de pesquisa da UNICAMP-Universidade Estadual
de Campinas. Nestas pesquisas são apresentados modelos matemáticos que simulam o
processo de degradação ambiental, que são estruturados através da equação de difusão-
advecção, nos quais buscam representar a dispersão de poluentes em diferentes ambientes
aquáticos regulares ou irregulares. Para o tratamento matemático utilizam do método das
diferenças finitas ou o método de elementos finitos. Neste sentido, traremos estudos, algumas
informações sobre a equação de difusão-advecção, com a intenção de construir um modelo
matemático teórico no domínio retangular. O tratamento numérico será realizado pelo
método das diferenças finitas para a discretização da equação, o método de Crank-Nicolson
para a parte das variáveis temporais, com objetivo de mostrar simulações computacionais
em ambiente Matlab da dispersão de poluentes em um ambiente aquático, baseados
no modelo matemático de Prestes (2011), fazendo apenas alterações nos parâmetro do
problema e preservando as mesmas condições de contorno, sendo elas, Robin Homogêneo e
Von Neumann.
AbstractThis paper presents some studies developed in the area of Biomathematics that has been
discussed by the research group of the UNICAMP-Universidade Estadual de Campinas.
These surveys are presented mathematical models that simulate the process of environ-
mental degradation, which are structured by diffusion-advection equation, in which seek to
represent the dispersion of pollutants in different regular or irregular aquatic environments.
For the mathematical treatment using the finite difference method or the finite element
method. In this sense, we will bring studies, some information on the diffusion-advection
equation, with the intention of building a theoretical mathematical model in rectangular
domain. The numerical treatment will be carried out by the finite difference method
for the discretization of the equation, the Crank-Nicolson method for the part of the
temporal variables, in order to show computer simulations in Matlab environment of
pollutant dispersion in an aquatic environment, based on the template mathematician
about citeonline, making only changes in the problem of parameter and preserving the
same boundary conditions, which were, Robin Homogeneous and Von Neumann.
Key-words: Diffusion, Advection, Dispersion, Pollutants, Matlab
Lista de ilustrações
Figura 1 – Gráfico de operadores de diferenças da série de Taylor. . . . . . . . . . 22
Figura 2 – John Crank. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
Figura 3 – Phyllis Nicolson. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
Figura 4 – O domínio retangular, da área que trabalhamos. . . . . . . . . . . . . . 25
Figura 5 – Sequência sucessiva de intervalos (malha). . . . . . . . . . . . . . . . . 26
Figura 6 – Pontos interiores da malha. . . . . . . . . . . . . . . . . . . . . . . . . 26
Figura 7 – Dispersão de poluente, sem influências do termo advectivo com u = 0 e
v = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
Figura 8 – Dispersão de poluentes, com influências do termo advectivo com u > 0
e v = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
Figura 9 – Dispersão de poluentes, com influências do termo advectivo com u < 0
e v = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Figura 10 – Dispersão de poluentes, com influências do termo advectivo com u = 0
e v > 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Figura 11 – Dispersão de poluentes, com influências do termo advectivo com u = 0
e v < 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
Figura 12 – Dispersão de poluente, com influências do termo advectivo com u = v. 34
Figura 13 – Dispersão de poluente, com influências do termo advectivo com u > v. 35
Figura 14 – Dispersão de poluente, com influências do termo advectivo com u < v. 35
Figura 15 – Dispersão de poluente, com influências do termo advectivo com u < 0 e
v > 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Figura 16 – Dispersão de poluente, com influências do termo advectivo com u > 0 e
v < 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Figura 17 – Dispersão de poluente, com influências do termo advectivo com u < v. 37
Figura 18 – Dispersão de poluente, com influências do termo advectivo com u > v. 37
Sumário
Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1 CONHECENDO A BIOMATEMÁTICA . . . . . . . . . . . . . . . . 13
1.1 Difusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2 Equação de Difusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.3 Métodos de Diferenças Finitas . . . . . . . . . . . . . . . . . . . . . 18
1.4 Série de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.5 Método de Crank-Nicolson . . . . . . . . . . . . . . . . . . . . . . . . 22
2 SIMULAÇÕES COMPUTACIONAIS E DISCUSSÕES . . . . . . . . 31
Conclusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
Referências . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
APÊNDICES 44
APÊNDICE A – CÓDIGO MATLAB . . . . . . . . . . . . . . . . . 45
ANEXOS 48
ANEXO A – CÓDIGO MATLAB USADO COMO MODELO DO
PRESTES . . . . . . . . . . . . . . . . . . . . . . . . . 49
Introdução
A terminologia "difusão"é utilizada em várias áreas do conhecimento: na sociologia,
economia e finanças, no crescimento populacional, opiniões, estimação de valores. Nas áreas
de exatas a equação de difusão é utilizada em pesquisas e estudos de áreas como física,
química, biologia, na engenharia mecânica, engenharia química, na engenharia elétrica e
outras.
De acordo com Pedron (2003) a difusão, de maneira geral, descreve o fenômeno de
um conjunto de partículas presentes em um sistema quando encaminha-se para o estado
de equilíbrio. Também pode ser descrita do ponto de vista físico como o “movimento
Browniano, definindo-se como movimento aleatório de partículas microscópicas imersas em
um fluido” (IFUSP, 2011, p. 2). Esse movimento irregular, de partículas macroscopicamente
pequenas (microscópicas) chamadas de moléculas ativas por Robert Brown, no entanto
maiores que as moléculas do fluido, está relacionada a fenômenos de transportes difusivos
em grande escala Michels (2014). Além disso, do ponto de vista da química, sua utilização
é de grande importância, quando trata-se no transporte de matéria através de propagação
no meio, devido ao movimento das partículas.
A dispersão das partículas em um dado sistema, é tido como difusão anômala,
sendo que a movimentação dessas partículas são aleatórias, permitindo que os fenômenos
de transporte sejam descritos através da equação de difusão.
A equação de difusão descreve o desenvolvimento de fenômenos físicos, biológicos,
químicos e outros que frequentemente ocorre na natureza. Ao observar esses fenômenos
podemos construir modelos matemáticos através de equações diferenciais parciais, também
aplicado na equação de difusão-advecção com o objetivo de representar esses fenômenos.
Dependendo da estrutura e construção do modelo para encontrar as soluções
dessa equação, podemos obter através das soluções analíticas ou encontrar por meio de
aproximações numéricas.
Assim, este trabalho foi estruturado de acordo com artigos, monografias, dis-
sertações e teses que se preocupam com a questão ambiental. Dentre esses trabalhos a
grande maioria vem sendo desenvolvido pelo grupo de pesquisa da UNICAMP. Estas
Introdução 12
pesquisas utilizam algumas ferramentas matemáticas, tais como o método das Diferenças
Finitas ou Elemento Finitos, a série de Taylor, o método de Crank-Nicolson e para as
aproximações numéricas usam o auxílio de programas computacionais que possibilita
apresentar simulações em ambiente Matlab de cenários que se degradam com o decorrer
do tempo, aplicam a equação de difusão-advecção para mostrar essa situação.
Por isso, nossa intenção foi desenvolver um estudo para conhecer estas pesquisas
relacionada a Biomatemática em especial, com a finalidade de construir de um modelo
teórico matemático que visa mostrar como a ocorre a dispersão de poluentes em um
ambiente aquático no domínio retangular.
No primeiro capitulo abordamos um pouco sobre o conceito de difusão o qual
usamos como referência Atkins e Jones (2006), Pedron (2003), em seguida a estrutura teórica
da equação de difusão-advecção posteriormente mostramos as ferramentas matemáticas
que escolhemos para o tratamento do modelo teórico, como o método das Diferenças
Finitas estruturada na série de Taylor e o método de Crank-Nicolson também usamos como
referência Ruggiero e Lopes (1996), Kaplan (1972), onde utilizamos como base os trabalhos
desenvolvidos pelo grupo de pesquisa da UNICAMP. No segundo capitulo apresentamos
algumas simulações computacionais do nosso modelo que usamos como apoio o trabalho
de Prestes (2011) para os estudos.
Afinal, apresentamos a conclusão, onde ressaltamos a importância deste trabalho
para minha formação.
1 Conhecendo a Biomatemática
O presente trabalho tem sua construção baseando-se em referências como artigos,
monografias, teses e dissertações como por exemplo Vanderlinde (2010), Melo (2011),
Pedron (2003), Israel (2011), Kirndges e Meyer (2011), Koga, Meyer e Tabares (2011),
Prestes (2011), Almeida (2010), Abreu (2009), Poletti (2009), Missio (2008), Wolmuth
(2009), onde alguns desses trabalhos realizados pelo grupo de pesquisa em Biomatemática
da UNICAMP – Universidade Estadual de Campinas. Este grupo desenvolve pesquisas
voltadas ao meio ambiente usando a modelagem matemática e simulações computacionais
em ambiente MATLAB, alguns desses estudos levam em consideração a influência do vento
no transportes da poluição, o efeito advectivo da equação de difusão-advecção, no domínio
bidimensional ou tridimensional sempre analisando cada caso, bem como a dispersão
de poluentes em meios aquáticos, que é um ramo de questão ambiental crescente em
discussões em nosso país. Esses modelos matemáticos visam principalmente representar
ações de recuperação dos cenários reais, através de cenários artificiais contribuindo para
um desenvolvimento sustentável do meio ambiente.
1.1 Difusão
As primeiras pesquisas experimentais, sobre a difusão foi realizada por Thomas
Graham, um químico escocês do século XIX que fez uma série de experiências sobre a
velocidade de efusão de gases. Descobriu-se também, que da mesma maneira, a difusão de
um gás varia aproximadamente com o inverso da raiz quadrada, Atkins e Jones (2006).
Décadas depois Adolf Fick (1855), usando as pesquisas de Graham, desenvolveu as leis
de difusão, observando a semelhança entre a difusão e a condução de calor de Fourier
em (1822) e a condução de eletricidade elétrica de Ohm em (1827) de acordo com Roque
(2014).
A difusão pode ser descrita pelo ponto de vista fenomenológico, referindo-se as
leis de Fick, definindo que as partículas transportam-se de áreas de maior concentração
para as áreas de menor concentração, que busca enfraquecer o gradiente de concentração
das partículas no ambiente Roque (2014).
Capítulo 1. Conhecendo a Biomatemática 14
Além disso, outra maneira descrever a difusão é abordado no ponto de vista
físico, atomística sobre o movimento Browniano e que foi mostrado bem antigamente, pela
primeira vez, com ênfase experimentalmente em:
Observe o que acontece quando os raios solares são admitidos em umprédio e lançam a luz sobre os lugares... Você verá uma infinidade depequenas partículas (grãos de poeira) se misturando em uma infinidadede maneiras... A dança destas partículas é uma indicação do movimentoda matéria que está escondida de nossa visão... Ele (o movimento) seoriginou com os próprios átomos [isto é, espontaneamente]. Então aquelespequenos corpos são colocados em movimento, removendo o ímpeto dosátomos. E, portanto, são postos em movimento pelo impacto de seusgolpes invisíveis (dos átomos nos grãos de poeira). Por sua vez, canhãocontra corpos um pouco maiores e assim por diante o movimento émontado a partir dos átomos e emerge gradualmente ao nível dos nossossentidos. Eu sei que os corpos estão em movimento. (Lucrécio, 60 a.capud IFUSP, 2011, p.2)
A descoberta do fenômeno Bowniano foi atribuída ao britânico escocês Robert
Brown (1773-1858), pois ele publicou o estudo em 1827. No entanto, segundo IFUSP
(2011), Brown não conseguiu descrever matematicamente esse fenômeno, só então em 1905
Albert Einstein explica a maneira certa de demostrar que o átomo existe.
Em suma, o nosso trabalho, utilizou os seguintes termos, sendo eles advecção
e dispersão. A advecção pode ser definida como transporte de matéria por um fluido,
onde pode haver movimentação do fluido. E a dispersão pode ser definida como a difusão
adicionada com dispersão hidrodinâmica, servindo no transporte de matéria no fluido de
forma mecânica de acordo com Dyminski (2006).
1.2 Equação de Difusão
A equação de difusão é formada por operações matemáticas definidas por equações
diferenciais parciais, a qual pode ser descrita por modelos matemáticos estruturados por
leis. Desse modo nos utilizamos as leis de Adolf Fick, onde usa o coeficiente de difusão,
para determinarmos a equação de difusão-advecção de acordo com Pedron (2003).
Observamos que alguns fenômenos de difusão segue a lei linear proposta por Fick,
assim temos:
J = −D∇ρ
• J representa a densidade de corrente;
Capítulo 1. Conhecendo a Biomatemática 15
• ∇ representa o vetor operador diferencial;
• D representa o coeficiente de difusão ou difusividade;
• ρ representa a densidade de elemento que se difundi.
Segue a demonstração de como determinar a equação de difusão, a partir das leis
de fick baseada na referência Pedron (2003).
Consideramos que a substância difundida não é nem absorvida nem emitida pelo
meio. Então a lei de conservação para esta substância implica uma equação de continuidade:
∂ρ
∂t+∇.J = 0
Combinando essas duas equações, chegamos à equação de difusão:
∂ρ
∂t= D∇2ρ
Se o sistema estiver sob a ação de uma força externa, ou arraste, a densidade de
corrente é:
J = −D∇ρ+ µFρ,
com µ representando a mobilidade. A equação de difusão com arraste é então
escrita como:
∂ρ
∂t= D∇2ρ− µ∇(Fρ)
A equação de difusão também é modificada se for possível a substância ser criada
(emitida) ou destruída (absorvida). Neste caso, a equação de continuidade é:
∂ρ
∂t+∇.J = δ,
na qual δ é a densidade da fonte, onde δ > 0 e δ < 0 estão associados à criação e
absorção de substância, respectivamente. Portanto, a correspondente equação de difusão,
não-homogênea, é dada por:
Capítulo 1. Conhecendo a Biomatemática 16
∂ρ
∂t= D∇2ρ+ δ.
Evidentemente, podemos ter os diversos fenômenos ocorrendo simultaneamente,
sendo descritos pela equação:
∂ρ
∂t= D∇2ρ− µ∇(Fρ) + δ.
Essa é a equação geral de difusão-advecção que é muito utilizada. Nos últimos
anos a equação de difusão vem contribuído em relação a transportes, onde usa a difusão
anômala em várias aplicações de modo geral, bem como, o transporte de fluidos em meios
porosos Ott et al. (1990), em análise batidas do coração em indivíduos saudáveis Peng et
al. (1993) estas são apenas dois exemplos mas ela pode ser usada em outras aplicações.
Atualmente há grande valorização sobre o desenvolvimento sustentável, por causa
da degradação do meio ambiente através de alguns fatores sendo eles, o crescente aumento
do desmatamento, a liberação de dióxido de carbono, a poluição do ar que é gerado
por queima de combustíveis fósseis como gasolina, diesel e carvão mineral, a poluição de
rios, lagos, mares e oceanos que são gerados por despejos de esgoto e acidente ambiental,
por exemplo, vazamento de petróleo, diesel e outros, como descreve o efeito estufa e
aquecimento global com isso tendo a possibilidade de alterações climáticas e entre outras
situações alterando o equilíbrio natural do ambiente.
Devido a esses fatores está sendo discutida em conferências por todo mundo e entre
países essas situações, uma reunião que ficou marcada e teve grande relevância nos últimos
anos foi a Conferência das Nações Unidas sobre o Meio Ambiente e Desenvolvimento
Sustentável(CNUDS) a mais conhecida como Rio+20 que aconteceu no Brasil em 2012,
foi realizada na cidade do Rio de Janeiro entre os dias 13 a 22 de junho, deixou para o
nosso país vários compromissos de desenvolvimento sustentável para os próximos anos, de
acordo com Brasil (2012).
Além dessas discussões, existe vários estudos voltados a preservar e diminuir esses
impactos, como é realizado pelo grupo de pesquisa em Biomatemática da UNICAMP, que
utilizam de ferramentas matemáticas como apresentadas anteriormente para descrever
esses processos principalmente em relação a degradação de ambientes aquáticos.
Capítulo 1. Conhecendo a Biomatemática 17
Por meio da construção do nosso trabalho, buscamos entender algumas ferramentas
matemáticas como o método de Diferenças Finitas, a Série de Taylor, o método de Crank-
Nicolson com o intuito de aplicá-las na equação de difusão-advecção, temos como objetivo
simular a dispersão de poluentes em um ambiente aquático, no caso de nosso modelo
teórico um domínio retangular.
Assim para apresentar esses procedimentos de dispersão de poluentes, usamos a
equação de difusão-advecção que foi aperfeiçoada através de uns modelos, tais como Okubo
e Levin (1980) para o coeficiente de difusão efetiva, como Edelstein-Keshet (1988) para
o decaimento, também conhecida como degradação o fenômeno chamado de fração das
partículas que sobre influências com meio externo de acordo com Prestes (2011), Marchuck
(1986) para o transporte advectivo u e v são vetores que são influenciados por fatores
externos principalmente ventos ou movimento aleatórios segundo Poletti (2009).
Indicamos a concentração de material impactante por C(x, y, t), sendo colocada
em determinada região num plano bidimensional (x, y, ) ∈ Ω ⊂ R2 verificada num instante
t ∈ (0, T ), assim reescrevemos a equação de difusão-advecção, como sendo:
∂C
∂t= div(α∇C)︸ ︷︷ ︸
Difusão
− div(CV)︸ ︷︷ ︸T ransporte
− σC︸︷︷︸Decaimento
+ f︸︷︷︸F onteP oluição
(1.1)
• C representa a concentração de material impactante;
• divα∇C representa a difusão efetiva;
• V = (u, v) representa as velocidades (advecção);
• σC representa a degradação (também conhecida como decaimento);
• f representa a fonte de poluição.
Essa é a equação que será utilizada para o desenvolvimento deste trabalho será
aplicado a métodos numéricos para a modelagem matemática e para as discretização das
equações, do domínio e para as variáveis temporais.
Capítulo 1. Conhecendo a Biomatemática 18
1.3 Métodos de Diferenças Finitas
Os problemas citados anteriormente, utilizam a equação de difusão-advecção, onde
muitas vezes caem equações diferenciais parciais não lineares que algumas apresentam
soluções e outras não apresentam uma solução analítica, o que necessita de um tratamento
diferenciado.
Assim percebemos que este trabalho vislumbra a simulação de um problema em
ambiental computacional, onde é utilizado o método das diferenças finitas que é conhecido
por simples aplicação na resolução de problemas de valor inicial ou de valor contorno,
podendo ser utilizado nas equações ordinárias ou em equações parciais.
O Método de Diferenças Finitas visa "transformar problemas formados por equa-
ções diferenciais em problemas envolvendo sistemas de equações algébricas, usando para
isto aproximações das derivadas que aparecem na equação, por diferenças finitas"conforme
Ruggiero e Lopes (1996). O método das diferenças finitas consiste na discretização de
derivadas, existentes na equação diferencial, onde essas derivadas são aproximadas pelas
diferenças finitas entre os valores da discretização.
De acordo com Ferreira; Ferreira, Kaibara e Navarro (2005) geralmente, a discre-
tização acontecesse em duas fases a primeira a discretização do domínio e a segunda a
discretização das equações.
Nesta primeira fase acontece a descrição numérica do domínio, contendo as posições
dos pontos da malha computacional onde buscamos encontrar a possível solução. Com isso,
consiste em dividir o domínio, formando vários subdomínio finitos, eles podem ser iguais
a uma malha uniforme, ou ser de tamanhos diferentes sendo uma malha não-uniforme.
Mas na maioria dos problemas usa-se a malha uniforme, em alguns casos, nos conceitos de
análise numérica possuem várias vantagens numéricas a utilização da malha não-uniforme.
Na segunda fase, é realizada a discretização das equações fornecendo as equações
de diferenças, com o objetivo de relacionar as variáveis dependentes nos pontos do domínio
discretizado.
Para construir as aproximações nas derivadas existentes na equação, do nosso
problema é necessário, usar os Operadores de Diferenças Finitas (ODF) que podem ser
encontrados pela expansão em Série de Taylor ou por interpolação polinomial. É possível
Capítulo 1. Conhecendo a Biomatemática 19
pela expansão em Série de Taylor obter os Operadores de Diferenças Atrasada, Operadores
de Diferenças Centrada e Operadores de Diferenças Adiantada.
Com isso os operadores de diferenças finitas faz a substituição da derivadas na
equação, dessa forma construindo um sistema de equações algébricas tornando os cálculos
simples para encontrar a solução e aplicá-lo na malha.
Neste trabalho usamos as aproximações nas derivadas de primeira e segunda
ordem. E a ferramenta matemática que utilizamos para encontrar essas aproximações das
derivadas será através da fórmula baseada na Série de Taylor.
1.4 Série de Taylor
Entende-se por série de Taylor, como sendo a expansão de séries de potências, que
possui uma função de derivadas que pode ser de qualquer ordem. De acordo com Kaplan
(1972). Seja f(h) a soma de uma série de potências cujo o intervalo de convergência é
a− r∗ < h < a+ r∗(r∗ > 0):
f(h) =∞∑
n=0Cn(h− a)n, a− r∗ < h < a+ r∗ (1.2)
Essa série denomina-se a série de Taylor de f(h) se os coeficientes Cn forem dados
pela regra:
C0 = f(a), C1 = f ′(a)1! , C2 = f ′′(a)
2! , Cn = f(n)n!
Temos então:
f(h) = f(a) + f ′(a)(h− a)1! + · · ·+ fn(n)(h− a)n
n! + · · ·
Kaplan (1972) apresenta o seguinte teorema, "toda série de potência com raio de
convergência não-nulo é a série de Taylor de sua soma". Demonstração: Seja f(h) dada por
Equação 1.2. Então por derivação pelo teorema descrito por Kaplan (1972), mostra que
pode-se derivar uma série de potência termo a termo dentro do intervalo de convergência,
se temos:
f(h) = C0 + C1(h− a) + C2(h− a)2 + · · ·+ Cn(h− a)n + · · ·
Capítulo 1. Conhecendo a Biomatemática 20
f ′(h) = C1 + 2C2(h− a) + · · ·+ nCn(h− a)n−1 + · · ·
f ′′(h) = 2C2 + 6C3(h− a) + · · ·+ n(n− 1)Cn(h− a)n−2 + · · · ...
fn(h) = n(n− 1) · · · 2Cn(n+ 1)n · · · 2Cn+1(h− a) + · · ·
Todas essas séries convergem em a − r∗ < h < a + r∗(r∗ > 0). Tomando agora
(h = a) obtemos:
f(a) = C0, f′(a) = C1, f
′′(a) = C2, · · · , fn(n) =n!Cn · · · ,
donde resulta que:
C0 = f(a) e Cn = fn(a)n! , n = 1, 2, · · · ,
Dessa forma, a série de Taylor pode ser escrita:
f(h) 6=+∞∑n=0
fn(a)n! (h− a)n = f(a) +f ′(a)(h−a) + f ′′(a)(h− a)2
2! + · · ·+ fn(a)n! (h−a) + · · ·
(1.3)
Com isso teremos as seguintes expressões, fazendo a substituição de (h = x−∆x)
e (h = x+ ∆x).
Para encontrar as aproximações das derivadas de primeira e segunda ordem,
a partir da série de Taylor utiliza-se o raio de convergência de ∆x, a função y será
infinitamente derivável como sendo variações de ”x” e a = x na equação 1.3.
Inicialmente iremos fazer a substituição e expandindo em y, substituindo (h =
x−∆x) temos:
f(h) = y(x−∆x) = y(x)+y′(x)(x−∆x−x)+y′′(x)2! (x−∆x−x)2−y
′′′(x)3! (x−∆x−x)3+· · ·+
(1.4)
Resolvendo:
y(x−∆x) = y(x)− y′(x)∆x+ y′′(x)2! ∆x2 − y′′′(x)
3! (∆x)3 + · · ·+ (1.5)
E de maneira análoga fazendo a substituição e expandindo em y, quando h =
x+ ∆x, teremos que:
f(h) = y(x+ ∆x) = y(x) + y′(x)(x+ ∆x− x) + y′′(x)2! (x−∆x− x)2 + y′′′(x)
3! (∆x)3 + · · ·+
Capítulo 1. Conhecendo a Biomatemática 21
Resolvendo:
y(x+ ∆x) = y(x) + y′(x)∆x+ y′′(x)2! ∆x2 − y′′′(x)
3! (∆x)3 · · ·+ (1.6)
Para encontrarmos a aproximação na ordem de ∆x3 para a primeira derivada
realizamos a subtração das equações 1.5 e 1.6, assim obtemos:
y(x−∆x) = y(x)− y′(x)∆x+ y′′(x)2! ∆x2 − y′′′(x)
3! (∆x)3 · · ·+
−
y(x+ ∆x) = y(x) + y′(x)∆x+ y′′(x)2! ∆x2 + y′′′(x)
3! (∆x)3 · · ·+
Resolvendo:
y(x−∆x)− y(x+ ∆x) = −2y′(x)∆x+ θ∆x3
⇒ y′(x) = y(x−∆x)− y(x+ ∆x)2∆x + θ∆x3
Da mesma maneira para determinar as aproximações na ordem de ∆x3 para a
segunda derivada, basta somar as equações 1.5 e 1.6, dessa forma termos que:
y(x−∆x) = y(x)− y′(x)∆x+ y′′(x)2! ∆x2 − y′′′(x)
3! + (∆x)3θ1∆x4 + · · ·+
+
y(x+ ∆x) = y(x) + y′(x)∆x+ y′′(x)2! ∆x2 + y′′′(x)
3! (∆x)3 + θ2∆x4 + · · ·+
Onde: θ3 = θ1 + θ2, resolvendo:
y(x−∆x) + y(x+ ∆x) = 2y(x) + 2y′′(x)2! + θ3∆x4
⇒ y(x−∆x) + y(x+ ∆x) = 2y(x) + y′′(x)∆x2 + θ3∆x4
⇒ −y(x−∆x)− 2y(x) + y(x+ ∆x)∆x2 θ3∆x4 = −y′′(x)(∆x2)
Capítulo 1. Conhecendo a Biomatemática 22
⇒ y(x−∆x)− 2y(x) + y(x+ ∆x)∆x2 + θ3∆x3∆x4 = y′′(x)
Dessa forma, podemos representar geometricamente a expansão da série de Taylor,
como determinamos nos cálculos acima os operadores da diferença atrasada e diferença
centrada, que também pode ser usada para determinar o operador de diferença adiantada,
como observamos na Figura 1.
Figura 1 – Gráfico de operadores de diferenças da série de Taylor.
1.5 Método de Crank-Nicolson
O método de Crank-Nicolson que utilizamos em nosso trabalho, foi desenvolvido
por dois físicos, sendo eles: John Crank e Phyllis Nicolson.
John Crank (1916-2006), recebeu o título de Bacharel e Mestre depois de Doutor
em 1953, posteriormente trabalhou na guerra balística como físico matemático em Cour-
taulds no Laboratório de Pesquisa Fundamental em (1945-1957), também foi professor
de matemática na Universidade de Brunel inicialmente na Brunel College, em Acton de
(1957-1981) Mathematics; Andrews (2000a).
Phyllis Nicolson é uma física que colaborou no desenvolvimento do método, o qual
trabalhamos. Estudou em Stockport High School recebendo o grau de Bacharel(1938),
Mestre (1939) e Doutorado em Física (1946) pela Universidade de Manchester e foi também
Capítulo 1. Conhecendo a Biomatemática 23
uma estudante de pesquisa (1945-1946) e bolsista de investigação (1946-1949) no Girton
College, Cambridge Mathematics; Andrews (2000a),
Figura 2 – John Crank.
fonte: Mathematics e Andrews (2000a)
Figura 3 – Phyllis Nicolson.
fonte: Mathematics e Andrews (2000b)
Enfim, ambos os pesquisadores John Crank e Phyllis Nicolson, tiveram em comum
grande reconhecimento de principal trabalho em analise numérica, em especial a solução
de problemas que envolve a equação do calor e as equações diferenciais parciais.
As fórmulas apresentadas anteriormente, da mesma maneira que usamos a expansão
em série de Taylor para obter os operadores de diferenças. Vamos determina a partir da
derivada de primeira ordem na variável t, buscando encontrar as aproximações para ∂C∂t,
com o erro na ordem de θ(∆t)2.
De acordo com Prestes (2011), aplicaremos as derivadas parciais em C(x, y, t) que
serão aproximadas em um mesmo instante, obtendo a discretização temporal em duas
etapas:
A primeira etapa será dividido um intervalo (0, T ] ⊂ R em nt subintervalos
[tn, tn+1] e ∆t = (T/nt).
Na segunda etapa, iremos considerar o método de Crank–Nicolson, que será
Capítulo 1. Conhecendo a Biomatemática 24
aplicado na equação diferencial através de derivadas parciais no ponto intermediário, na ∂C∂t
nos pontos t(n+∆t/2), e representamos por C(xi, yi, tn) sendo Cni e C(xi, yi, tn+∆t/2) sendo
C(n+1/2)i .
Para desenvolver o método, utilizaremos os cálculos auxiliares abaixo, são eles:
C(n)i = C
(n+1/2)i − ∆t
2∂C
∂t+ (−∆t)2
2∂C2
∂t2+ θ1(−∆t)3 (1.7)
C(n+1)i = C
(n+1/2)i + ∆t
2∂C
∂t+ (∆t)2
2∂C2
∂t2+ θ2(∆t)3 (1.8)
Realizando a operação de subtração das equações 1.7 e da 1.8:
C(n)i = C
(n+1/2)i − ∆t
2∂C
∂t+ (−∆t)2
2∂C2
∂t2+ θ1(−∆t)3
−
C(n+1)i = C
(n+1/2)i + ∆t
2∂C
∂t+ (∆t)2
2∂C2
∂t2+ θ2(∆t)33
Resolvendo:
C(n)i − C(n+1)
i = −∆t∂C∂t
+ θ3(∆t)2
⇒ ∂C
∂t= C
(n+1)i − C(n)
i
∆t θ3(∆t)2
Somando as equações 1.7 e da 1.8, no instante t = t(n+1/2), temos que:
C(n)i = C
(n+1/2)i − ∆t
2∂C
∂t+ (−∆t)2
2∂C2
∂t2+ θ2(−∆t)3
+
C(n+1)i = C
(n+1/2)i + ∆t
2∂C
∂t+ (∆t)2
2∂C2
∂t2+ θ2(∆t)3
Resolvendo:
C(n)i + C
(n+1)i = 2C(n+1/2)
i + θ3(∆t)2
⇒ C(n+1/2)i = C
(n)i + C
(n+1)i
2 + θ3(∆t)2
Com estas equações realizamos a discretização temporal, utilizando as derivadas
parciais em relação x e y que serão aproximadas, com os resultados encontrados, ou
seja, trabalhamos com um método estável e vantajoso para aproximações numéricas na
utilização computacional de acordo com Carnahan, Luther e Wilkes (1969) apud Prestes
(2011).
Capítulo 1. Conhecendo a Biomatemática 25
O modelo matemático que apresentamos estará utilizando a equação de difusão-
advecção apresentada anteriormente, o qual estaremos discretizando pelo método das
diferenças finitas e Crank-Nicolson na parte temporal.
A região do domínio (Ω), será uma área retangular de comprimento c e largura r,
em relação as bordas usamos as condições de contorno baseadas no trabalho de Prestes
(2011), sendo elas Robin Homogêneo e Von Neumann.
Figura 4 – O domínio retangular, da área que trabalhamos.
Onde será, divido o eixo das abscissas em subintervalos [0, c] em nx vários subin-
tervalos, assim ∆x = cnx, da mesma forma dividimos o eixo das ordenadas [0, r] em ny
subintervalos, temos que ∆y = rny.
Construindo a malha na região r, considerando ny a quantidade de divisores de
"r"e nx a quantidade de divisores em "c", temos:
Capítulo 1. Conhecendo a Biomatemática 26
Figura 5 – Sequência sucessiva de intervalos (malha).
Considerando as condições iniciais do problema iremos manter as bordas como
referência de Prestes (2011), logo a quantidade de pontos na malha será de nny = ny + 1
pontos por colunas e nnx = nx + 1 pontos por linha, totalizando (tp = (nny ∗ nnx)).
Apresentado um ponto arbitrário Ci, mostrarmos a segui na malha os pontos interiores de
um determinado domínio, Figura 6, que em nosso caso é uma situação retangular:
Figura 6 – Pontos interiores da malha.
Determinados os pontos interiores na malha, utilizaremos a Equação 1.1, reescrevendo-
a de forma que obteremos as derivadas de segunda ordem em relação a x e a y, tendo que
o V vetor é constante, onde por u é a direção em relação ao eixo x e v é a direção ao eixo
Capítulo 1. Conhecendo a Biomatemática 27
y, assim teremos que a Equação 1.1 ficará:
∂C
∂t= α
(∂2c
∂x2 + ∂2C
∂y2
)︸ ︷︷ ︸
Difusão
−u(∂C
∂x
)− v
(∂C
∂y
)︸ ︷︷ ︸
T ransporte
− σC︸︷︷︸Decaimento
+ f︸︷︷︸F onteP luição
(1.9)
Com isso, desenvolvemos o método de Crank–Nicolson, usando os cálculos auxilia-
res, e as equações parciais de primeira e segunda ordem que demostramos anteriormente
pela série de Taylor, aplicado na malha de nós em f(xi, yi, t(n+∆t/2)), é indicado por
f (n+1/2) para encontrar a solução aproximada, de acordo com os cálculos das equações
vistas anteriormente estabelecemos:
∂C
∂t(ti, t(n+∆/2)) = C
(n+1)i − C(n)
i
∆t∂2C
∂x2 (xi, t(n+∆t/2)) ∼=C
(n+1/2)i+nny − 2C(n+1/2)
i + C(n+1/2)i−nny
∆x2 ∗
∂2C
∂y2 (yi, t(n+∆/2)) ∼=C
(n+1/2)i+1 − 2C(n+1/2)
i + C(n+1/2)i−1
∆y2
∂C
∂x(xi, t(n+∆/2)) ∼=
C(n+1/2)i+nny − C
(n+1/2)i−nny
2∆x ∗
∂C
∂y(yi, t(n+∆/2)) ∼=
C(n+1/2)i+1 − C(n+1/2)
i−12∆y
Observação: Temos que (∗) são as derivadas parciais, em relação a ”x” e de acordo
com a Figura 6. O elemento anterior ”ai” será i− nny e posterior i+ nny
C(n+1)i − Cn
i
∆t∼=
α
C(n+1/2)i+nny − 2C(n+1/2)
i + C(n+1/2)i−nny
∆x2
+C(n+1/2)
i+1 − 2C(n+1/2)i + C
(n+1/2)i−1
∆y2
−u
C(n+1/2)i+nny − C
(n+1/2)i−nny
2∆x
− vC(n+1/2)
i+1 − C(n+1/2)i−1
2∆y
− σC(n+1/2)i + f (n+1/2)
Usando a equação acima, vamos trabalhar com a Equação 1.9 na aplicação do
Método de Crank-Nicolson, fazendo as devidas substituições:
C(n+1)i − C(n)
i
∆t∼= α
C
(n+1)i+nny + C
(n)i+nny
2 − 2C(n+1)i + 2C(n)
i
2 +C
(n+1)i−nny + C
(n)i−nny
2∆x2
+
Capítulo 1. Conhecendo a Biomatemática 28
α
C
(n+1)i+1 + C
(n)i+1
2 − 2C(n+1)i + 2C(n)
i
2 + C(n+1)i−1 + C
(n)i−1
2∆y2
−
u
(C(n+nny)1+nny + C
(n)i+1)− (C(n+1)
i−nny + C(n)i−nny)
4∆x
− v(C(n+1)
i+1 + C(n+1)i+1 )− (C(n+1)
i−1 + C(n)i−1)
4∆y
−
σ
(C(n+1)i + C
(n)i )
2
+ fn+1/2
Na equação a seguir, à variável ∆t irá multiplicar os termos do lado direito da
equação, na mesma oportunidade vamos resolver algumas frações dos termos na equação,
dessa forma obtemos:
C(n+1)i − C(n)
i∼= α∆t
C(n+1)i+nny + C
(n)i+nny − (2C(n+1)
i + 2C(n)i ) + (C(n+1)
i−nny + C(n)i−nny)
2∆x2
+
α∆tC(n+1)
i+1 + C(n)i+1 − (2C(n+1)
i + 2C(n)i ) + (C(n+1)
i−1 + C(n)i−1)
2∆y2
+
u∆t(C(n+1)
i+nny + C(n)i+nny)− (C(n+1)
i−nny + C(n)i−nny)
4∆x
−
v∆t(C(n+1)
i+1 + C(n)i+1)− (C(n+1)
i−1 + C(n)i−1)
4∆y
−
σ∆tC(n+1)
i + C(n)i
2
+
∆tf (n+1/2)
Com isso, agrupamos os termos semelhantes da equação acima em relação a n+ 1
e n, assim:
C(n+1)i −C(n)
i∼= α∆t
C(n+1)i+nny − 2C(n+1)
i + C(n+1)i−nny
2∆x2
+α∆tC(n+1)
i+1 − 2C(n+1)i + C
(n+1)i−1
2∆y2
+
α∆tC(n)
i+nny − 2C(n)i + C
(n)i−nny
2∆x2
+ α∆tC(n)
i+1 − 2C(n)i + C
(n)i−1
2∆y2
−
Capítulo 1. Conhecendo a Biomatemática 29
u∆tC(n+1)
i+nny − C(n+1)i−nny
4∆x
− v∆tC(n+1)
i+1 − C(n+1)i−1
4∆y
−
u∆tC(n)
i+nny − C(n)i−nny
4∆x
− v∆tC(n)
n+1 − C(n)i−1
4∆y
−
σ∆tC(n+1)
i
2
− σ∆tC(n)
i
2
+
∆tf (n+1/2)
A equação anterior fizemos manipulações algébricas, assim a equação ficará:
C(n+1)i − α∆t
C(n+1)i+nny − 2C(n+1)
i + C(n+1)i−nny
2∆x2
α∆tC(n+1)
i+1 − 2C(n+1)i + C
(n+1)i−1
2∆y2
+
u∆tC(n+1)
i+nny − C(n+1)i−nny
4∆x
+ v∆tC(n+1)
i+1 − C(n+1)i−1
4∆y
+ σ∆tC(n+1)
i
2
=
C(n)i + α∆t
C(n)n+nny − 2C(n)
i + C(n)i−nny
2∆x2
α∆tC(n)
i+1 − 2C(n)i + C
(n)i−1
2∆y2
−
u∆tC(n)
i+nny − C(n)i−nny
4∆x
− v∆tC(n)
i+1 − C(n)i−1
4∆y
−
σ∆tC(n)
i
2
+ f (n+1/2)∆t
Dessa forma, para resolver a equação acima simplificamos os termos semelhantes,
com o intuito de deixar em evidência os termos em relação a n+ 1 e n, com isso temos:
(− α∆t
2∆x2 −u∆t4∆x
)C
(n+1)i−nny +
(− α∆t
2∆x2 + v∆t4∆y
)C
(n+1)i−1 +
(1 + α∆t
∆x2 + α∆t∆y2 + σ∆t
2
)C
(n+1)i +
(− α∆t
2∆y2 + v∆t4∆y
)C
(n+1)i+1 +
(− α∆t
2∆x2 + u∆t4∆x
)C
(n+1)i+nny =
(α∆t2∆x2 + u∆t
4∆x
)C
(n)i−nny +
(α∆t2∆y2 + v∆t
4∆y
)C
(n)i−1 +
(1− α∆t
∆x2 −α∆t∆y2 −
σ∆t2
)C
(n)i +
Capítulo 1. Conhecendo a Biomatemática 30
(α∆t2∆y2 −
v∆t4∆y
)C
(n)i+1 +
(α∆t2∆x2 −
u∆t4∆x
)C
(n)i+nny+
f (n+1/2)∆t
A partir da equação acima podemos estruturar um sistema linear formado por
equações e incógnitas relacionando certos nós de uma malha, de acordo com Prestes (2011),
de uma região regular ou irregular.
Respeitando a malha, Prestes (2011) construiu um programa em ambiente Matlab,
com o intuito de representar como a poluição se dispersa em uma região por meio
de simulações computacionais. Para realizar essas representações através de simulações
computacionais, deve-se estabelecer algumas condições de contorno tais como Dirichlet,
Neumann e Robin podendo aplicar nas bordas das regiões estudas.
Para apresentar as simulações nesde trabalho, não analisaremos as condições de
contorno do problema, simplesmente estaremos assumindo as mesmas condições apresenta-
das por Prestes (2011).
2 Simulações Computacionais e Discussões
Como o presente trabalho foi baseado em outros estudos que trazem contribuições
importantes para a questão ambiental, um deles usamos como modelo matemático a
referencia de Prestes (2011), assim como também o programa em ambiente Matlab que
utilizou em sua dissertação.
Realizamos simulações no programa apresentado por Prestes (2011) em ambiente
Matlab, fazendo modificações nos parâmetros do programa, utilizamos as mesmas condições
de contorno, sendo elas, Robin homogêneo e Von Neumann para apresentar como uma
fonte poluição se dispersa em uma região retangular.
A fonte de poluição que utilizamos em nossos estudos, foi empregada por exemplo
um cano de esgoto que estourou sobre uma região, dessa forma a poluição é libera de
forma continua sobre um o ambiente aquático.
De acordo com a equação de difusão-advecção com região Ω, apresentaremos assim
os gráficos abaixo que representam o processo de dispersão de poluentes, com 168 iterações,
sendo uma fonte poluente pontual localizada no centro da região estuda, consideramos essa
região retangular como por exemplo um plano cartesiano assim determinarmos a direção
da dispersão de poluentes.
Primeiramente considerando o vetor advectivo ~V = (u, v) como elemento nulo.
Nestas condições, a poluição se dispersão se espalha, sem influência de velocidades ou
qualquer corrente.
Capítulo 2. Simulações Computacionais e Discussões 32
Figura 7 – Dispersão de poluente, sem influências do termo advectivo com u = 0 e v = 0.
Observamos que o movimento apresentado na acima não é circular, isto ocorrem
devido as condições de contorno apresentado por Prestes (2011).
As figuras 8 e 9 a seguir, iremos considerar os vetores u > 0 e v = 0, depois
u < 0 e v = 0, assim percebermos que quando u representa a valores maiores ou menores
que zero e v = 0, a poluição sofre influências de vento, ou fluxo de água entre outras
influências, ocorrendo dispersão direcionada no sentido relacionado ao eixo x para direita
e para esquerda.
Figura 8 – Dispersão de poluentes, com influências do termo advectivo com u > 0 e v = 0.
Capítulo 2. Simulações Computacionais e Discussões 33
Figura 9 – Dispersão de poluentes, com influências do termo advectivo com u < 0 e v = 0.
O processo de dispersão das figuras 8 e 9 deveria seguir para direita ou para
esquerda tal mudança se dá pelas condições de contorno.
Situação semelhante acontece nas figuras 10 e 11 seguintes onde u = 0 e v > 0,
logo após u = 0 e v < 0 nessas condições onde v assume valores maiores ou menores que
zero e u = 0 a poluição se dispersa, com influências externas do meio onde se localiza a
poluição como por exemplo correntezas, em relação ao eixo y positivo e negativo.
Figura 10 – Dispersão de poluentes, com influências do termo advectivo com u = 0 e
v > 0.
Capítulo 2. Simulações Computacionais e Discussões 34
Figura 11 – Dispersão de poluentes, com influências do termo advectivo com u = 0 e
v < 0.
Nas figuras 12, 13 e 14 seguintes podemos observar que os valores de u e v são
negativos, porém na Figura 12 u e v são valores iguais, também negativos. Podemos
concluir que sendo ambos valores diferentes negativos, nas três figuras a seguir ocorre se
dispersão de poluentes com influências externas do meio para o mesmo sentido, podendo
ser por exemplos chuvas, observamos o que modifica são as velocidades, pois a cada figura
diminui a dispersão de poluente de forma direcionada para o eixo y positivo:
Figura 12 – Dispersão de poluente, com influências do termo advectivo com u = v.
Capítulo 2. Simulações Computacionais e Discussões 35
Figura 13 – Dispersão de poluente, com influências do termo advectivo com u > v.
Figura 14 – Dispersão de poluente, com influências do termo advectivo com u < v.
Nas próximas figuras15 e 16 temos que u > 0 e v < 0, e em situação contraria que
seria u < 0 e v > 0, assim verificamos a dispersão de poluente ocorre em sentidos opostos
em relação ao eixo y, mostrando como a poluição se espalha com a influência do meio
externo, como ventos, chuva e entre outros.
Capítulo 2. Simulações Computacionais e Discussões 36
Figura 15 – Dispersão de poluente, com influências do termo advectivo com u < 0 e v > 0.
Figura 16 – Dispersão de poluente, com influências do termo advectivo com u > 0 e v < 0.
As últimas duas figuras 17 e 18 possuem a mesma características da figura anterior
pois u e v são maiores que zero e positivos, porém possuem valores diferentes, semelhante
as figuras anteriores, elas mostram a dispersão do poluente para o mesmo sentido, no
entanto a cada figura a dispersão aumenta em relação ao eixo y. Com influências do meio,
como por exemplo ventos ou movimento ondulatórios, observe:
Capítulo 2. Simulações Computacionais e Discussões 37
Figura 17 – Dispersão de poluente, com influências do termo advectivo com u < v.
Figura 18 – Dispersão de poluente, com influências do termo advectivo com u > v.
Observamos em todas as figuras 2 até 18, quando v 6= 0 o direcionamento contrário
ao plano coordenado em relação aos eixos x e y, assumindo u como eixo x e v como eixo
y. Notamos alteração na direção de v, pois quando v > 0 a dispersão direciona para
a parte negativa do gráfico e quando v < 0 direção é para a parte positiva do gráfico,
não descobrimos se tal situação ocorre por questão do processo ou se por erro na parte
computacional.
Enfim, essas simulações computacionais da equação de difusão-advecção pode ser
Capítulo 2. Simulações Computacionais e Discussões 38
utilizada para mostrar, tipos fenômenos que trabalhamos relacionado a Biomatemática e
entres outros fenômenos aplicados em outras áreas da ciência.
Conclusão
Considerando um dos objetivos que tentamos alcançar para o presente estudo
foi construir e fazer a programação de um modelo matemático para desenvolver algumas
simulações computacionais de como a poluição se dispersa em uma região retangular.
Assim definimos estudar um modelo matemático de Prestes (2011) com o intuito de realizar
algumas simulações computacionais fazendo modificações nos parâmetros do problema
para conhecer como podemos utilizar a equação de difusão-advecção para representar
processos físico na propagação de poluentes em ambientes aquáticos. Para isto, com base
em Prestes (2011) que utilizou com o software Matlab, mostrando o desenvolvimento da
dispersão de poluentes na superfície de uma região retangular.
Nosso trabalho objetivou em termos, um primeiro contato com estudos que
analisam o comportamento da concentração de poluente na superfície de uma região
aquática, sobre a influência de fatores externos, tais como, vento, chuvas, movimento
ondulatórios, fluxo das águas até correntezas entre outras influências, para isto, utilizamos
o termo advectivo no domínio retangular Ω para trabalhar com a equação de difusão-
advecção. Neste trabalho destacamos a superfície interna os pontos interiores do domínio
Ω, e não as bordas do domínio que foi estudado, portanto resolvemos deixar para fazer os
estudos do tratamento das bordas por meio das condições de contorno, tais como: Diriclet,
Von Neumann e Robin para futuros trabalhos, como a elaboração de um artigo ou futura
dissertação de mestrado.
Este trabalho proporcionou a possibilidade de pesquisas e estudos matemáticos
em Análise Numérica para conhecer o processo de dispersão de poluentes em ambientes
áquaticos, usando o método de Diferenças Finitas, a série de Taylor e o método de Crank-
Nicolson, com o intuito de verificar o comportamento advectivo presente na equação de
difusão-advecção, com isto pudemos conhecer os conteúdos matemáticos que são pouco
abordados no curso de Licenciatura em Matemática, com isso, ampliando o campo de
estudos para a continuação e processos de formação na área de matemática, se possível
em um programa de mestrado.
Como o objetivo do presente trabalho foi conhecer os estudos e pesquisas na
Conclusão 40
área de Biomatemática, em especial do grupo da UNICAMP que abordam o método
das diferenças finitas e como ocorre a dispersão de poluentes em ambientes áquaticos
em domínios regulares ou irregulares, que apresentam simulações computacionais em
ambiente Matlab, como é o desenvolvimento do método e o resultado apresentados pelos
pesquisadores estudados, como queríamos demostrar concluímos que o nosso objetivo
principal foi alcançado.
41
Referências
ABREU, L. C. Influência de Poluentes sobre Macroalgas na Baía de Sepetiba, RJ:Modelagem Matemática, Análise Numérica e Simulações Computacionais. Dissertação(Mestrado) — Universidade Estadual de Campinas, Instituto de Matemática, Estatística eComputação Científica, São Paulo, 2009.
ALMEIDA, J. R. F. Modelagem Matemática e Simulação Computacional para análise deDispersão de Poluente em um trecho do Rio Paraíba do Sul. Dissertação (Mestrado) —Universidade Estadual de Campinas, Instituto de Matemática, Estatística e ComputaçãoCientífica, São Paulo, 2010.
ATKINS, P.; JONES, L. Princípios de Química: Questionando a Vida Moderna e o MeioAmbiente. 3. ed. [S.l.]: Bookman, 2006.
BRASIL. Rio+20. [S.l.], 2012. Disponível em: <http://www.rio20.gov.br/sobre_a_rio_mais_20.html>. Acesso em: 03 nov de 2013.
CARNAHAN, B.; LUTHER, H. A.; WILKES, J. O. Applied Numerical Methods,. NewYork: John Wiley & Sons, Inc., 1969.
DYMINSKI, A. S. Contaminação de Solos e Águas Subterrâneas. [S.l.], 2006. Disponívelem: <definir>. Acesso em: 03 nov de 2013.
EDELSTEIN-KESHET, L. Mathematical Models in Biology. 1. ed. México: McGraw-Hill,Inc., 1988.
FERREIRA, V. G.; KAIBARA, M. K.; NAVARRO, H. A. Modelagem Matemática eSimulação Númerica em Dinâmica dos Fluidos. São Paulo: SBMAC, 2005.
IFUSP. Movimento Browniano. São Paulo, 2011. Disponível em: <http://web.if.usp.br/ifusp/files/Browniano-L.pdf>. Acesso em: 03 de nov de 2014.
ISRAEL, V. P. Análise Bayesiana de Processos de Difusão Multivariados com ObservaçõesDiscretas. Dissertação (Mestrado) — Universidade Federal do Rio de Janeiro Instituto deMatemática Departamento de Métodos Estatísticos, 2011.
KAPLAN, W. Cálculo Avançado. [S.l.]: Edgard Blücher, 1972.
KIRNDGES, A.; MEYER, F. da C. A. Modelagem e simulação computacional dadispersão de poluentes no lago de manso-mt, com uso da equação de difusão-advecção.Biomatemática, São Paulo, v. 21, p. 193–208, 2011.
KOGA, M. T.; MEYER, F. da C. A.; TABARES, P. C. C. Dinâmica populacionalinterativa da mosca-dos-chifres (haematobia irritans) na presença de um predador:Simulações computacionais. Biomatemática, São Paulo, v. 21, p. 71–86, 2011.
MARCHUCK, I. G. Mathematicals models in environmental problems, studies inmathematics and its applications. Elsevier Science Publishers, v. 16, p. 217, 1986.
Referências 42
MATHEMATICS, S. of; ANDREWS, S. Statistics University of S. John Crank. [S.l.], 2000.Disponível em: <http://www-history.mcs.st-andrews.ac.uk/Biographies/Crank.htm>.
MATHEMATICS, S. of; ANDREWS, S. Statistics University of S. Phyllis Nicolson. [S.l.],2000. Disponível em: <http://www-history.mcs.st-andrews.ac.uk/Biographies/Nicolson.html>.
MELO, K. J. M. de. Aplicação do Método das Diferenças Finitas Explícito na Solução daEquação do Calor para o Caso Transiente e Unidimensional. Dissertação (Mestrado) —Universidade Federal Rural do Semi-Árido Campus Angicos Curso de Ciência e Tecnologia,2011.
MICHELS, F. S. Aplicaçõeses da Equacão de Difusão: Vínculo Geomêtrico e DifusãoIônica. Dissertação (Mestrado) — Universidade Estadual de Maringá, 2014.
MISSIO, M. Modelos de EDP integrados à Lógica Fuzzy e Métodos Probalisticos notratamento de incertezas: uma aplicação à febre aftosa em bovinos. Tese (Doutorado) —Universidade Estadual de Campinas, Instituto de Matemática, Estatística e ComputaçãoCientífica, São Paulo, 2008.
OKUBO, A.; LEVIN, S. A. Diffusion and Ecological Problems: Mathematical Models. 2.ed. Berlim: Springer, 1980.
OTT, A. et al. Anomalous diffusion in "living polymers": a genuine levy flight? Phys, p.65:2202, 1990.
PEDRON, I. T. Estudos em Difusão Anômala. Dissertação (Mestrado) — UniversidadeEstadual de Maringá Departamento de Física, 2003.
PENG, C. K. et al. Long range anticorrelations and non-gaussian behavior of theheartbeat. Phys, p. 70:1343, 1993.
POLETTI, E. C. C. Dispersão de poluente em Sistema de Reservatório: ModelagemMatemática e Simulação Computacional utilizando-se Aproximação Numérica e ConjuntosFuzzy. Tese (Doutorado) — Universidade Estadual de Campinas, Faculdade de EngenhariaElétrica e de Computação, São Paulo, 2009.
PRESTES, M. F. B. Dispersão de material impactante em meio aquático: ModeloMatemático, aproximação numérica e simulação computacional- Lagoa do TaquaralCampinas, SP. Dissertação (Mestrado) — Universidade Estadual de Campinas, Institutode Matemática, Estatística e Computação Científica, São Paulo, 2011.
ROQUE, A. Biofisica II - Difusão. [S.l.], 2014. Disponível em: <http://sisne.org/Disciplinas/Grad/Biofisica2FisMed/aula1.pdf>. Acesso em: 03 nov de 2013.
RUGGIERO, M. A. G.; LOPES, V. L. R. Cálculo Numérico: Aspactos Teóricos eComputacionais. 2. ed. São Paulo: Person Makron Books, 1996.
VANDERLINDE, J. B. Aplicação do Método das Diferenças Finitas na Resoluções deEquações Diferenciais Ordinárias. Dissertação (Mestrado) — Universidade de Estado deMato Grosso - Campus Universitário de Sinop, 2010.
Referências 43
WOLMUTH, L. D. Modelagem e simulações do comportamento evolutivo de poluentesem corpos aquáticos de grande extensão: o caso da represa do rio Manso. Dissertação(Mestrado) — Universidade Estadual de Campinas Instituto de Matemática, Estatística eComputação Científica, São Paulo, 2009.
Apêndices
APÊNDICE A – Código Matlab
1 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % Cenario : Dominio retangular %
3 % Condicoes de Fronteira : Von Neumann e Robin homogeneo nas bordas %
4 % Direcao do Vento: Basta alterar os%
5 % componentes dos ventos : Nos parametros do problema %
6 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 clear all
8 % entrada de parametros do problema
9 alfa =0.7e -3; sigma =0.2e -4;
10 u= -0.09e -2; v= -0.19e -2; f=1;
11 % parametros do dominio
12 l=2.0; h=1.0; tf =600;
13 % parametros da discretizacao
14 nx =40; dx=l/nx; ny =40; dy=h/ny;
15 nnx=nx +1; nny=ny +1; nn=nnx*nny;
16 nt =168; dt=tf/nt;
17 % Ncleo de Peclet
18 ux=u*dx/alfa;vy=v*dy/alfa;
19 % condicao inicial
20 czero=zeros(nn ,1);
21 % calculos auxiliares
22 ddx=dx*dx; ddy=dy*dy;
23 adtx=alfa*dt/ddx; adty=alfa*dt/ddy;
24 sdt2=sigma*dt /2;
25 udt=u*dt /(4* dx); vdt=v*dt /(4* dy);
26 %
27 dpe =1+ adtx+adty+sdt2;
28 dpd =1-adtx -adty -sdt2;
29 dsipe=-adty /2- vdt; dsipd=adty /2+ vdt;
30 dsspe=-adty /2+ vdt; dsspd=adty /2- vdt;
31 dside=-adtx /2- udt; dsidd=adtx /2+ udt;
32 dssde=-adtx /2+ udt; dssdd=adtx /2- udt;
33 % preparacao das matrizes e do termo - fontes
34 me= sparse (nn);md= sparse (nn);
35 b=zeros(nn ,1);
APÊNDICE A. Código Matlab 46
36 % preenchimento de me , me e b
37 for i=1: nn
38 me(i,i)=dpe; md(i,i)=dpd;
39 end
40 for i=1:nn -1
41 me(i,i+1)=dsspe; md(i,i+1)=dsspd;
42 me(i+1,i)=dsipe; md(i+1,i)=dsipd;
43 end
44 for i=1:nn -nny
45 me(i,i+nny)=dssde; md(i,i+nny)=dssdd;
46 me(i+nny ,i)=dside; md(i+nny ,i)=dsidd;
47 end
48 % correces devidas s bordas
49 % borda inferior
50 me (1 ,2)=-adty; md (1 ,2)=adty;
51 for j=1: nnx -1
52 in=j*nny +1;
53 me(in ,in -1) =0; md(in ,in -1) =0;
54 me(in ,in +1)=-adty; md(in ,in +1)=adty;
55 end
56 % borda superior
57 for j=1: nnx -1
58 in=j*nny;
59 me(in ,in -1)=-adty; md(in ,in -1)=adty;
60 me(in ,in +1) =0; md(in ,in +1) =0;
61 end
62 me(nn ,nn -1)=-adty; md(nn ,nn -1)=adty;
63 % borda esquerda
64 for i=1: nny
65 me(i,i+nny)=-adtx; md(i,i+nny)=adtx;
66 end
67 % borda direita
68 for i=1: nny
69 in=(nnx -1)*nny+i;
70 me(in ,in -nny)=-adtx; md(in ,in -nny)=adtx;
71 end
72 % mudanas devidas condicao de Robin
73 % ndices dos ns do vertedouro
74 for k=1:4
APÊNDICE A. Código Matlab 47
75 in=nn -(k+2);
76 me(in ,in -nny)=-adtx /2;
77 md(in ,in -nny)=adtx /2;
78 me(in ,in)=1+ adtx+adty+sdt2+dt/dx -u*dt /(2* alfa);
79 md(in ,in)=1-adtx -adty -sdt2 -dt/dx+u*dt /(2* alfa);
80 end
81 % termo independente
82 in =19* nny +22; b(in)=f*dt;
83 % resoluo repetida do sistema
84 ver=zeros(nny ,nnx);
85 for it =1: nt
86 c=me\( md*czero+b);
87 czero=c;
88 % para poder ver
89 for i=1: nny
90 for j=1: nnx
91 ind =(j -1)*nny+i;
92 ver(nny +1-i,j)=c(ind);
93 end
94 end
95 contour (ver ,30) ,grid
96 if (mod(it ,150))==0 pause
97 colorbar ;
98 end
99 pause (0.1)
100 end
Anexos
ANEXO A – Código Matlab usado como
modelo do Prestes
1 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % Cenrio : lagoa retangular %
3 % Condies de Fronteira : Von Neumann e Robin homogneo no vertedouro %
4 % Direo do Vento: Noroeste , para o Sudeste , basta alterar os%
5 % componentes dos ventos : Sudeste , nos parmetros do problema %
6 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 clear all
8 % entrada de parametros do problema
9 alfa =0.125e -3; sigma =0.25e -8;
10 u= -0.19e -2; v=0.09e -2; f=1;
11 % parametros do dominio
12 l=2.0; h=1.0; tf =200;
13 % parametros da discretizacao
14 nx =32; dx=l/nx; ny =20; dy=h/ny;
15 nnx=nx +1; nny=ny +1; nn=nnx*nny;
16 nt =750; dt=tf/nt;
17 % Ncleo de Peclet
18 ux=u*dx/alfa;vy=v*dy/alfa;
19 % condicao inicial
20 czero=zeros(nn ,1);
21 % calculos auxiliares
22 ddx=dx*dx; ddy=dy*dy;
23 adtx=alfa*dt/ddx; adty=alfa*dt/ddy;
24 sdt2=sigma*dt /2;
25 udt=u*dt /(4* dx); vdt=v*dt /(4* dy);
26 %
27 dpe =1+ adtx+adty+sdt2;
28 dpd =1-adtx -adty -sdt2;
29 dsipe=-adty /2- vdt; dsipd=adty /2+ vdt;
30 dsspe=-adty /2+ vdt; dsspd=adty /2- vdt;
31 dside=-adtx /2- udt; dsidd=adtx /2+ udt;
32 dssde=-adtx /2+ udt; dssdd=adtx /2- udt;
33 % preparacao das matrizes e do termo - fontes
ANEXO A. Código Matlab usado como modelo do Prestes 50
34 me= sparse (nn);md= sparse (nn);
35 b=zeros(nn ,1);
36 % preenchimento de me , me e b
37 for i=1: nn
38 me(i,i)=dpe; md(i,i)=dpd;
39 end
40 for i=1:nn -1
41 me(i,i+1)=dsspe; md(i,i+1)=dsspd;
42 me(i+1,i)=dsipe; md(i+1,i)=dsipd;
43 end
44 for i=1:nn -nny
45 me(i,i+nny)=dssde; md(i,i+nny)=dssdd;
46 me(i+nny ,i)=dside; md(i+nny ,i)=dsidd;
47 end
48 % correces devidas s bordas
49 % borda inferior
50 56
51 me (1 ,2)=-adty; md (1 ,2)=adty;
52 for j=1: nnx -1
53 in=j*nny +1;
54 me(in ,in -1) =0; md(in ,in -1) =0;
55 me(in ,in +1)=-adty; md(in ,in +1)=adty;
56 end
57 % borda superior
58 for j=1: nnx -1
59 in=j*nny;
60 me(in ,in -1)=-adty; md(in ,in -1)=adty;
61 me(in ,in +1) =0; md(in ,in +1) =0;
62 end
63 me(nn ,nn -1)=-adty; md(nn ,nn -1)=adty;
64 % borda esquerda
65 for i=1: nny
66 me(i,i+nny)=-adtx; md(i,i+nny)=adtx;
67 end
68 % borda direita
69 for i=1: nny
70 in=(nnx -1)*nny+i;
71 me(in ,in -nny)=-adtx; md(in ,in -nny)=adtx;
72 end
ANEXO A. Código Matlab usado como modelo do Prestes 51
73 % mudanas devidas condicao de Robin
74 % ndices dos ns do vertedouro
75 for k=1:4
76 in=nn -(k+2);
77 me(in ,in -nny)=-adtx /2;
78 md(in ,in -nny)=adtx /2;
79 me(in ,in)=1+ adtx+adty+sdt2+dt/dx -u*dt /(2* alfa);
80 md(in ,in)=1-adtx -adty -sdt2 -dt/dx+u*dt /(2* alfa);
81 end
82 % termo independente
83 in =5* nny +10; b(in)=f*dt;
84 % resoluo repetida do sistema
85 ver=zeros(nny ,nnx);
86 for it =1: nt
87 c=me\( md*czero+b);
88 czero=c;
89 % para poder ver
90 for i=1: nny
91 for j=1: nnx
92 ind =(j -1)*nny+i;
93 ver(nny +1-i,j)=c(ind);
94 end
95 end
96 contour (ver ,30) ,grid
97 if (mod(it ,150))==0 pause
98 end
99 pause (0.1)
100 end