Embed Size (px)
Citation preview
UNIVERSIDADE FEDERAL DO CEARÁ
CENTRO DE TECNOLOGIA
DEPARTAMENTO DE ENGENHARIA HIDRÁULICA E AMBIENTAL
MESTRADO EM ENGENHARIA CIVIL
ÁREA DE CONCENTRAÇÃO EM RECURSOS HÍDRICOS
FRANCISCO DAS CHAGAS AZEVEDO DOS REIS
MODELAGEM MATEMÁTICA E COMPUTACIONAL DA CONTAMINAÇÃO DE
AQUÍFEROS COM USO DE MÉTODOS NUMÉRICOS SEM MALHA
FORTALEZA – CEARÁ 2014
ii
FRANCISCO DAS CHAGAS AZEVEDO DOS REIS
MODELAGEM MATEMÁTICA E COMPUTACIONAL DA CONTAMINAÇÃO DE
AQUÍFEROS COM USO DE MÉTODOS NUMÉRICOS SEM MALHA
Dissertação submetida à Coordenação do Curso de Pós-Graduação em Engenharia Civil, com área de concentração em
Recursos Hídricos, como parte dos requisitos para a obtenção do título de Mestre em Engenharia Civil.
Orientador: Prof. PhD. Marco Aurélio
Holanda de Castro
FORTALEZA – CEARÁ
2014
iii
Dados Internacionais de Catalogação na Publicação
Universidade Federal do Ceará
Biblioteca de Pós-Graduação em Engenharia - BPGE
R31m Reis, Francisco das Chagas Azevedo dos. Modelagem matemática e computacional da contaminação de aquíferos com uso de métodos
numéricos sem malha / Francisco das Chagas Azevedo dos Reis. – 2014.
117 f. : il. color., enc. ; 30 cm.
Dissertação (mestrado) – Universidade Federal do Ceará, Centro de Tecnologia, Departamento
de Engenharia Hidráulica e Ambiental, Programa de Pós-Graduação em Engenharia Civil, Fortaleza,
2014.
Área de Concentração: Recursos Hídricos
Orientação: Prof. Dr. Marco Aurélio Holanda de Castro.
1. Recursos hídricos. 2. Equações diferenciais parciais. 3. Difusão. I. Título.
CDD 627
iv
FRANCISCO DAS CHAGAS AZEVEDO DOS REIS
MODELAGEM MATEMÁTICA E COMPUTACIONAL DA CONTAMINAÇÃO DE
AQUÍFEROS COM USO DE MÉTODOS NUMÉRICOS SEM MALHA
Dissertação submetida à Coordenação do Curso de Pós-Graduação em Engenharia
Civil, com área de concentração em Recursos Hídricos, como parte dos requisitos para a obtenção do título de
Mestre em Engenharia Civil. Área de concentração: Recursos Hídricos.
v
À Deus. A minha mãe, esposa, filhas e
familiares em geral, os quais servem de
base para tudo que realizo.
vi
AGRADECIMENTOS
Agradeço em especial à minha esposa Patrícia Nunes dos Reis, pelo seu
apoio irrestrito na confecção desse trabalho, as minhas filhas que toleraram as
minhas ausências. Ao meu orientador professor Marco Aurélio Holanda de Castro,
que sempre foi incansável e paciente, atento ao rigor, sempre nos ajudou, em todos
os momentos, na confecção desse trabalho. Agradeço também aos colegas João
Marcelo Costa Barbosa, José Valmir Farias Maia Junior e Nelci Rones Pereira de
Sousa por nos ajudar, com o SCILAB.
Ao Instituto Federal do Piauí (IFPI), Universidade Federal do Piauí (UFPI),
Universidade Estadual do Piauí (UESPI) e a Universidade Federal do Ceará (UFC),
pela celebração do convênio, no qual deu origem ao MINTER UFC/IFPI. A todos os
professores do Departamento de Engenharia Hidráulica e Ambiental (DEHA), em
especial, aos professores que participaram do MINTER: Marco Aurélio Holanda de
Castro, Francisco de Assis de Souza Filho, Iran E. Lima Neto, José Capelo Neto,
Ticiana M. Carvalho Studart, José Nilson B. Campos e John Kenedy de Araújo pelo
inestimável aprendizado.
vii
RESUMO
Em muitos problemas da natureza e em uma diversidade enorme de áreas do conhecimento, existe a necessidade real de modelarmos fenômenos existentes. Em Ciências como Matemática, Física, Química, Biologia, Economia e nas Engenharias,
de uma maneira geral, é comum por parte dos pesquisadores, o uso de modelos e simulações, às quais, quase sempre, envolvem taxas, princípios e leis, regidos por Equações Diferenciais. Problemas envolvendo movimento de fluidos, intensidade de
corrente elétrica, propagação de calor, crescimento populacional, entre muitos outros, são exemplos clássicos de aplicações de modelos regidos por Equações Diferencias, às quais, podem ser diferenciadas quanto ao tipo em Equações Diferenciais
Ordinárias (EDO) e Equações Diferenciais Parciais (EDP). Nas primeiras, a função a ser determinada depende de uma única variável independente, enquanto nas segundas, ocorre a dependência de duas ou mais variáveis independentes. Acontece
é que em uma grande variedade de problemas da natureza, as equações não possuem soluções bem comportadas, analíticas e, dessa maneira, faz-se necessário o conhecimento de métodos numéricos, tais como, Diferenças Finitas, Elementos
Finitos, Elementos de Contorno, entre outros, os quais necessitam da discretização do domínio e, portanto da criação de uma malha (MESH), com fórmulas interativas
para se estimar uma solução e minimizar o erro da aproximação. Nesse sentido, a
proposta desse trabalho é utilizar um método numérico bastante eficaz e independente de malha, denominado método sem malhas (MESHLESS), mas
especificamente o método de Kansas, o qual lança mão de Funções de Base Radial (Radial Basis Functions – RBF), ou simetria radial, da distância entre um ponto central
do domínio da função e um ponto genérico do domínio. A função interpoladora de base radial, também depende de um parâmetro de forma “c” a ser encontrado. Mas a questão preponderante é: como determinar um parâmetro de forma “c” ótimo, que
possa oferecer uma solução consistente, reduzindo o resíduo e, portanto o erro existente? Para tanto, modelou-se um problema de contaminação de aquífero fazendo uso da equação de difusão, comparando o resultado de sua solução
analítica, com a solução numérica obtida através do método numérico sem malhas e com o parâmetro de forma simulado e otimizado por meio da plataforma SCILAB (versão 5.4.1).
Palavras-Chave: Métodos Sem Malhas, Modelagem Matemática, Contaminação de Aquíferos, Equação da Difusão, Equações Diferenciais.
viii
ABSTRACT
In many problems of nature and a huge diversity of knowledge areas, there is a real
need we model existing phenomena. Sciences like Mathematics, Physics, Chemistry,
Biology, Economics and in Engineering, in general, is common among the
researchers, the use of models and simulations, which almost always involve fees,
principles and laws, governed by Differential Equations. Problems involving fluid
motion, intensity of electric current, heat propagation, population growth, among many
others, are classic examples of applications of models governed by Differential
Equations, which can be differentiated as to type in Ordinary Differential Equations
(ODE) and Partial Differential Equations (PDE). In the first, the function to be
determined depends on a single variable, while in the second, the dependence of two
or more independent variables occurs. Happens is that in a wide variety of problems of
nature, the equations do not have well-behaved, analytic and thus solutions, it is
necessary the knowledge of numerical methods such as Finite Differences, Finite
Elements, Boundary Elements, among others, which require the discretization of the
domain and therefore the creation of a mesh (MESH), with interactive formulas for
estimating a solution and minimize the error of approximation. In this sense, the
purpose of this work is to use a very efficient and independent of mesh numerical
method, called method without mesh (MESHLESS), but specifically the method of
Kansas, which makes use of Radial Basis Function (Radial Basis Functions - RBF) or
radial symmetry, the distance between central point of the domain of the function and
a generic point of the domain. The interpolating radial basis function also depends on
a shape parameter "c" to be found. But the overriding question is how to determine a
shape parameter "c" great, we can provide a consistent solution, reducing waste and
therefore the existing error? For both, modeled itself a problem of contamination of the
aquifer by making use of the diffusion equation, comparing the results of its analytical
solution with the numerical solution obtained by numerical method without mesh and
parameter simulated shape and optimized by SCILAB platform (version 5.4.1).
Keywords: Meshless, Mathematical Modeling, Contamination of Aquifers, Diffusion
Equation, Differential Equations.
ix
LISTA DE ILUSTRAÇÕES
Figura 2.1 – Ciclo Hidrológico ........................................................................................ 3
Figura 2.2 – Raio de convergência de uma série ....................................................... 16
Figura 2.3 – Barra de propagação de calor................................................................. 22
Figura 4.1 – Parâmetro de forma utilizado "c=0,755929" , pontos do domínio, com
número de pontos N=5................................................................................................. 45
Figura 4.2 – Parâmetro de forma utilizado "c=0,755929" , pontos do domínio, com
número de pontos N=7. .............................................................................................. 46
Figura 4.3 – Parâmetro de forma utilizado "c=0,755929" , pontos do domínio, com
número de pontos N=10. ............................................................................................ 46
Figura 4.4 – Definição das condições de contorno da EDO....................................... 50
Figura 4.5 – Discretização do domínio da EDO. ........................................................ 50
Figura 4.6 – Definição dos coeficientes da EDO. ....................................................... 51
Figura 4.7 – Dados de entrada para a otimização do parâmetro de forma c. ........... 51
Figura 4.8 – Gráfico de resíduo do contorno e determinação de c = 68.21 como sendo
o c - ótimo, via método multiquadrático direto. ........................................................... 52
Figura 4.9 – Gráfico de resíduo do domínio e determinação de c=15.81, como sendo o
c - ótimo, via método multiquadrático direto................................................................ 52
Figura 4.10 – Comparação entre solução numérica e solução analítica, pontos do
domínio, com c=15.81, otimizado pelo método multiquadrático direto. ..................... 53
Figura 4.11 – Definição das condições de contorno da EDO, para uso do método
inverso. ......................................................................................................................... 53
Figura 4.12 – Gráfico de resíduo e determinação de c=70.61, como sendo o c - ótimo
do contorno, via método multiquadrático inverso. ...................................................... 54
Figura 4.13 – Gráfico de resíduo e determinação de c=16.81, como sendo o c - ótimo
do domínio, via método multiquadrático inverso. ........................................................ 54
Figura 4.14 – Comparação entre solução numérica e analítica, pelo método inverso
c=16.81. ........................................................................................................................ 55
Figura 4.15 – Solução analítica da Equação da Difusão Pura para os instantes t=0s;
t=0.05s; t=0.10s com N=50 pontos, x (Km). ................................................................ 62
Figura 4.16 – Comparação solução analítica com solução numérica da Equação da
Difusão Pura, pontos do domínio, para os instantes t=0s; t=0.05s; t=0.10s com N=5
pontos, x (Km). ............................................................................................................. 63
x
Figura 4.17 – Comparação solução analítica com solução numérica da Equação da
Difusão Pura, pontos do domínio, para os instantes t=0s; t=0.05s; t=0.10s com N=20
pontos, x (Km). ............................................................................................................. 63
Figura 4.18 – Comparação solução analítica com solução numérica da Equação da
Difusão Pura, pontos do domínio, c-ótimo 0.3, para os instantes t=0s; t=0.05s; t=0.10s
com N=40 pontos. ........................................................................................................ 64
Figura 4.19 – Gráfico comparativo de c e do somatório do quadrado do resíduo .... 64
Figura 4.20 – Fluxograma simplificado da otimização do c. ....................................... 73
Figura 4.21 – Gráfico comparativo da solução numérica e da solução analítica da
Equação de Difusão-Advecção com t=15 min e N=100 pontos. ................................ 75
Figura 4.22 – Gráfico comparativo da solução numérica e da solução analítica da
Equação de Difusão-Advecção com t=15 min e N=200 pontos. ................................ 75
Figura 4.23 – Gráfico comparativo entre a solução numérica e a solução analítica da
Equação de Difusão-Advecção, pontos do domínio, nos instantes t=12min, t=24min,
t=36min. O parâmetro de forma, c = 0.149254, foi otimizado pelo SCILAB 5.4.1..... 76
Figura 4.24 – Gráfico comparativo do parâmetro de forma e do resíduo da equação de
Difusão-Advecção.. ...................................................................................................... 76
Figura A.1 – Expansão do vetor u ............................................................................. 82
Figura A.2 – Mudança de sentido do vetor u ............................................................ 83
Figura A.3 – Gráfico da função delta de Dirac de . .................................... 86
Figura A.4 – Gráfico da função delta de Dirac de , quando ............. 86
Figura B.1– Sobre o SCILAB 5.4.1 .............................................................................. 87
Figura B.2 – Ambiente, plataforma SCILAB 5.4.1 ...................................................... 88
Figura B.3 – Acessando o aplicativo Scinotes do SCILAB 5.4.1 ............................... 88
Figura B.4 – Aplicativo SciNotes do SCILAB 5.4.1 ..................................................... 89
Figura B.5 – Gerando uma matriz 3x3 no SCILAB 5.4.1 ............................................ 89
Figura B.6 – Determinando a dimensão de uma matriz A, no SCILAB 5.4.1 ............ 90
Figura B.7 – Cálculo da transposta de uma matriz A, no SCILAB 5.4.1................... 91
Figura B.8 – Cálculo da inversa de uma matriz A, no SCILAB 5.4.1 ......................... 92
Figura B.9 – Esquema de uso do comando “for” no SCILAB 5.4.1 ........................... 93
Figura B.10 – Aplicação do comando “for” no SCILAB 5.4.1 ..................................... 94
Figura B.11 – Aplicação do comando “while” no SCILAB 5.4.1 ................................. 95
Figura B.12 – Navegador de ajuda no SCILAB 5.4.1, possui comandos específicos
para os mais variados ramos da Matemática. ............................................................ 96
xi
LISTA DE TABELAS
Tabela 2.1 - Água disponível na terra ........................................................................... 4
Tabela 2.4.1– Valores de Difusividade Térmica ......................................................... 25
Tabela 3.1.1– Funções de base radial ........................................................................ 27
Tabela 4.1.1– Discretização dos pontos do domínio da EDO de 1ª ordem............... 40
Tabela 4.1.2– Matriz R das distâncias ........................................................................ 40
Tabela 4.1.3– Matriz FI = ϕ, função de interpolação, ponto a ponto. ......................... 41
Tabela 4.1.4– Matriz DFI das derivadas da função FI, ponto a ponto. ...................... 41
Tabela 4.1.5 – Matriz LFI, soma das matrizes FI e DFI .............................................. 42
Tabela 4.1.6– Matriz de colocação, A, é o produto da inversa de λ pela matriz F .... 42
Tabela 4.1.7– Inversa da matriz de colocação A ........................................................ 43
Tabela 4.1.8– Resumo dos dados da solução da EDO de 1ª ordem via método sem
malhas .......................................................................................................................... 43
Tabela 4.1.9– Componentes do vetor resíduo, VR, da EDO de 1ª ordem ................ 44
Tabela B.1 - Operadores básicos do SCILAB ............................................................. 94
xii
LISTA DE ABREVIATURAS E SIGLAS
EDO Equação diferencial ordinária
EDP Equação diferencial parcial
MESHLESS Método Livre de Malhas
DBO Demanda Bioquímica de oxigênio
OD Oxigênio Dissolvido
MVP Método da Variação de Parâmetros
LI Linearmente independente
LD Linearmente dependente
SCILAB Software matemático
RBF Função de base radial
VR Vetor resíduo
MQ Multiquadrática direta
MQI Multiquadrática inversa
xiii
LISTA DE SÍMBOLOS
L Demanda bioquímica de oxigênio
Constante de desoxigenação do efluente
Fator de integração
Solução particular
Função complementar
Raio de convergência
Constante de difusividade térmica
A Matriz de colocação
Função de interpolação
e Matrizes de armazenamento dos pontos de discretização do domínio
R Matriz das distâncias
c Parâmetro de forma
Vetor que contém os coeficientes da solução aproximada
Número de difusão
E Coeficiente de dispersão
C(x,t) Concentração
Massa total das partículas normais à área da seção transversal
Número de Courant
Velocidade de escoamento
xiv
SUMÁRIO
1. INTRODUÇÃO ........................................................................................................ 1
2. REVISÃO BIBLIOGRÁFICA .................................................................................. 2
2.1 Ciclo Hidrológico ............................................................................................... 2
2.2 Distribuição de água no planeta ....................................................................... 4
2.3 Contaminação de aquíferos ............................................................................. 5
2.3.1 Propagação de poluentes no meio aquático............................................. 6
2.4 Equações diferenciais....................................................................................... 9
2.4.1 Classificação de equações diferenciais quanto ao tipo ............................ 9
2.4.2 Classificação de equações diferenciais quanto à ordem ......................... 9
2.4.3 Classificação de equações diferenciais quanto à linearidade................ 10
2.4.4 Solução de equação diferencial linear de 1ª ordem ............................... 10
2.4.5 Equações diferenciais de 2ª ordem ......................................................... 11
2.4.6 Equações diferenciais de ordem superior ............................................... 13
2.4.7 Aproximação em série de potências para uma equação diferencial linear
de segunda ordem ................................................................................................. 14
2.4.8 Série de potências ................................................................................... 15
2.4.9 Série de Fourier ....................................................................................... 20
2.4.10 Equações diferenciais parciais e o método de separação de variáveis. 21
3. METODOLOGIA.................................................................................................... 26
3.1 O método sem malhas: Funções de Base Radial ......................................... 27
3.2 Operadores lineares aplicados em uma EDP................................................ 28
3.3 Método de Kansa ............................................................................................ 29
3.4 Equação geral apresentada no SCILAB ........................................................ 31
4. SIMULAÇÕES E RESULTADOS ......................................................................... 34
4.1 Aplicação do método sem malhas (MESHLESS) em uma EDO via Método de
Kansa. ....................................................................................................................... 34
4.1.1 Vetor Resíduo - VR .................................................................................. 44
4.2 Aplicação do método sem malhas (MESHLESS) em uma EDO de 2ª ordem
com coeficientes constantes via método de Kansa................................................. 47
4.2.1 Fluxograma do processo de uso do método sem malhas...................... 49
4.3 Equação da difusão pura................................................................................ 55
4.4 Equação Diferencial Parcial da Difusão Advecção ....................................... 65
4.4.1 Problema envolvendo a equação da difusão advecção ......................... 66
xv
4.4.2 Abordagem Kansa na equação de difusão advecção ............................ 67
4.4.3 Matrizes da equação de difusão advecção ............................................. 70
4.4.4 Estabilidade do sistema Número de Courant ......................................... 74
5. CONCLUSÕES: .................................................................................................... 77
REFERÊNCIAS BIBLIOGRÁFICAS ........................................................................... 79
APÊNDICE A – NOCÕES DE DEPENDÊNCIA LINEAR E INDEPENDÊNCIA LINEAR
E FUNÇÃO DE IMPULSO DE DIRAC. .................................................................... 82
A.1 – Conjunto linearmente dependente (LD) ......................................................... 82
A.2 – Conjunto linearmente independente (LI) ........................................................ 83
A.3 – Aplicação Linear (Transformação Linear) ...................................................... 84
A.4 – Aplicação linear do vetor nulo entre espaços vetoriais. ................................ 84
A.5 – Funçõesde Impulso ......................................................................................... 85
APÊNDICE B – Breve exposição sobre o SCILAB ................................................ 87
B.1 – INTRODUÇÃO ................................................................................................ 87
B 1.1 – Criação de Matrizes .................................................................................. 89
B 1.2 – Dimensão de Matrizes .............................................................................. 90
B 1.3 – Transposição de Matrizes ........................................................................ 91
B 1.4 Inversa de Matrizes...................................................................................... 92
B 1.5 – Comandos para Iterações (for e while) .................................................... 93
B 1.6 – Comandos e funções pré-definidos no SCILAB ...................................... 96
APÊNDICE C - ROTINA DO PROBLEMA 4.1 DA EDO DE 1ª ORDEM. ................. 97
APÊNDICE D- ROTINA DO PROBLEMA DA EDO DE 2ª ORDEM NÃO
HOMOGÊNEA. ........................................................................................................... 100
APÊNDICE E- ROTINA DO PROBLEMA DA EQUAÇÃO DE DIFUSÃO-PURA. . 108
APÊNDICE G- ROTINA DA APROXIMAÇÃO DA SOLUÇÃO NUMÉRICA COM A
SOLUÇÃO ANALÍTICA DO PROBLEMA DA EQUAÇÃO DE DIFUSÃO-ADVECÇÃO
CONSIDERANDO O TEMPO DE 15 min. ................................................................ 115
1
1. INTRODUÇÃO
A água é um recurso natural de suma importância para todos os seres
vivos. Dessa forma, é importante nós cuidarmos bem de tal elemento tão precioso.
Nesse sentindo, é fundamental que tenhamos alguns conhecimentos básicos acerca
dessa substância e dos elementos que podem degradá-la tornando-a imprópria para
o consumo.
Há algum tempo a sociedade de maneira geral vem se preocupando com
a questão da contaminação de aquíferos pela ação antrópica, crescimento
desordenados das cidades, poluição de mananciais, entre outros problemas. Com
isso, vários estudiosos de diversas áreas, entre as quais, os de Engenharia de
Recursos Hídricos, têm se debruçados com intento de descobrir soluções para tais
casos.
A modelagem matemática e computacional tem servido para nortear uma
imensa gama de problemas, simulando situações, nas quais, estudos in locu, sejam
dispendiosos e muitas vezes inacessíveis. Nesse sentindo, várias técnicas de
modelagem, ferramentas computacionais e softwares, são utilizadas no meio
científico para facilitar e descrever modelos matemáticos, haja vista grande parte
desses fenômenos, envolverem Equações Diferenciais Ordinárias (EDO) e
Equações Diferenciais Parciais (EDP). Em tais equações, nem sempre é fácil
descobrir uma solução analítica e, em problemas de Engenharia muitas vezes, nem
existe uma solução desse tipo, para isso é notório a importância dos Métodos
Numéricos, os quais possam estimar uma solução “confiável”, mas que não
inviabilizem o problema com erros devido a processos de simplificações ou devido
ao truncamento e discretização de dados feitos pelo computador.
Entre os vários métodos capazes de solucionar Equações Diferencias,
destacamos o Método das Diferenças Finitas, o de Elementos Finitos, Elementos de
Contorno e o Método Livre de Malhas (MESHLESS). Nossa proposta é fazer uma
comparação entre uma solução analítica da Equação de Dispersão de um poluente,
ao longo de um trecho de um aquífero, determinada através da aplicação de série de
Fourier e o Método Livre de Malhas, mais especificamente as funções de base
radial, e o método de Kansas.
2
2. REVISÃO BIBLIOGRÁFICA
A ciência que estuda a água é denominada HIDROLOGIA, do grego
HYDOR que significa “água” e LOGOS, que é traduzido como “ciência”. Nesse
sentido, define-se Hidrologia como:
Hidrologia é a ciência que trata da água na Terra, sua ocorrência, circulação
e distribuição, suas propriedades físicas e químicas, e sua reação com o
meio ambiente, incluindo sua relação com as formas vivas relacionada com
toda a água da Terra, sua ocorrência, distribuição e circulação, suas
propriedades físicas e químicas, seu efeito sobre o meio ambiente e sobre
todas as formas da vida. (CHOW, 1959, p.1).
2.1 Ciclo Hidrológico
Ciclo Hidrológico é o nome dado à circulação da água em nosso planeta,
desde o momento de seu estado natural, até a sua evaporação, subida a atmosfera
em forma de vapor, e descida seguindo o caminho inverso, em forma de
precipitação, na qual ao tocar o solo, pode escoar superficialmente ao encontro de
oceanos, lagos, lagoas e rios, ou, ainda, infiltrar-se no solo e percolar até camadas
mais profundas e ser armazenadas como água subterrânea, podendo a posteriori
fluir como manancial, para que todo o processo ocorra novamente.
A força necessária para a movimentação e continuação do ciclo
hidrológico, ad eternum, é a energia solar, a qual aquece o ar, o solo e a superfície
da água, além das folhas das plantas ocasionando o processo conhecido por
evapotranspiração, o qual de maneira simplificada, podemos dizer que é a
evaporação das gotículas de água contidas nas folhas dos vegetais. Dessa forma,
esse aquecimento por parte do Sol também é causador das massas de ar, às quais
fazem circular as nuvens e provocam chuvas.
3
O ciclo hidrológico é um fenômeno complexo composto de vários
processos: precipitação, interceptação, transpiração, evaporação, infiltração,
percolação, armazenamento e escoamento (MIJARES, 1992, p.17).
Fonte: USGS
Figura 2.1 – Ciclo Hidrológico
4
2.2 Distribuição de água no planeta
Apesar do ciclo hidrológico, ser um fenômeno complexo e
permanentemente movido pela força solar, a distribuição de água no planeta é
bastante irregular. Os fatores dessa irregularidade são diversos, entre os quais
temos a própria variedade de clima no planeta, relevo, solo, cobertura vegetal, entre
outros fatores.
A água é a substância mais abundante na biosfera terrestre e a mais
importante, junto com o oxigênio para manutenção da vida no planeta, contudo a
maior parte dela está nos oceanos e mares, os quais correspondem acerca de 97%
da água do planeta. Dos três porcentos restante, grande parte em torno de 97% está
inacessível, em camadas profundas do subsolo acerca de 800 metros, em aquíferos
subterrâneos; em forma de gelo nas camadas polares e em regiões de altas
montanhas, tem-se 2,4%; algo como 0,6% está na atmosfera e apenas o número de
0,3%, em termos de volume, têm-se quatro milhões de quilômetros cúbicos, os quais
podem ser aproveitados de alguma forma pelos homens (PHILIPPI, et. al. 2004).
Este último, por sua vez, é constantemente degradado pela ação antrópica, pelo uso
de defensivos agrícolas, esgotos domésticos e industriais, entre outros efluentes não
tratados.
Tabela 2.1 - Água disponível na terra
Fonte: PHILIPPI, et. al.(2004)
LOCAL VOLUME (x103km3) PORCENTAGEM DA ÁGUA TOTAL
Lagos de água doce 125 0,009
Rios 1,25 0,000000919
Umidade do solo 65 0,000047
Água subterrânea 8.250 0,0060
Lagos-salinos e mares interiores
105 0,0007
Atmosfera 13 0,000009
Calotas de gelo polares, geleiras e neve
29.200 0,02147
Oceano e mares 1.320.000 0,9706
TOTAL 1.360.000 100
5
2.3 Contaminação de aquíferos
Muito embora a água seja um recurso fundamental para a subsistência
de todos os seres vivos, a poluição de aquíferos, rios, lagos e mananciais têm
atingido níveis alarmantes, os quais despertaram, nos últimos anos, a preocupação
da sociedade de modo geral e, em particular do meio acadêmico, onde observamos
um número muito grande de cientistas e estudantes que se debruçam e tentam
através de suas pesquisas, dar respostas no sentindo de mitigar a deletéria dos
recursos hídricos disponíveis.
Nesse sentindo, é importante entendermos os processos de
contaminação e de poluição. O termo poluição vem do latim polluere, que significa
sujar, manchar, e dessa forma, vem sendo empregado indiscriminadamente, sem o
devido cuidado. Haja vista a ideia no sentindo estritamente estético que uma água
turva, ou com determinados microrganismos, seja uma água poluída, ou imprópria
para o uso de qualquer natureza é errônea em sua essência. Para aclarar os fatos, é
importante entendermos a distinção entre poluição e contaminação, a saber: Nessa
perspectiva:
Entende-se por poluição da água a alteração de suas características por
quaisquer ações ou interferências, sejam elas naturais ou provocadas pelo
homem. Essas alterações podem produzir impactos, estéticos, fisiológicos
ou ecológicos (BRAGA et. al. 2005, p. 82).
Podemos dizer que a contaminação é a transferência de elementos
nocivos à saúde. Ou seja, o fato de a água estar contaminada, por elementos
patogênicos à saúde humana, não implica, necessariamente, em poluição e morte
daquela biosfera local, assim como, o fato da água está com alterações em suas
características devidas a processos, como calor excessivo, ou a turbidez devido ao
carreamento de sedimentos não são, indubitavelmente, dolosos à saúde humana.
No tocante a aquíferos, sabe-se que apesar dos mesmos de uma maneira
geral possuírem um processo natural de autodepuração, através de bactérias e
fungos que se alimentam de poluentes orgânicos biodegradáveis. Estes, por sua
vez, em quantidade excessiva e pela ação das próprias bactérias podem ocasionar a
morte dos peixes, devido ao aumento da Demanda Bioquímica de Oxigênio (DBO) e,
6
por conseguinte, a redução do Oxigênio Dissolvido (OD). Além disso, existem
poluentes orgânicos do tipo recalcitrantes ou refratários, os quais não são
biodegradáveis, ou tem uma taxa de degradação muito lenta, entre eles podemos
citar os defensivos agrícolas, detergentes sintéticos e o petróleo. Ademais
compostos inorgânicos como os metais pesados, a exemplo do mercúrio usado em
atividades de garimpo e do chumbo cujos potenciais carcinogênicos são altamente
comprometedores.
Para aprofundarmos nas características de aquíferos, ressaltamos a
existência de dois tipos de aquíferos subterrâneos, os artesianos (ou confinados) e
os freáticos. Ambos mesmo dotados de sua proteção natural, no subsolo, não estão
isentos de sofrerem com os efeitos de vários tipos de contaminação, as mais
comuns, são de agrotóxico, combustível, fossas, cemitérios e à intrusão salina, para
os costeiros.
Enfim, existe uma série de efluentes danosos à saúde da população os
quais poderiam ser evitados com o uso de práticas mais conscientes de preservação
dos recursos hídricos e de preservação de nascentes.
Nesse sentindo, a modelagem computacional-matemática, surge como
ferramenta importante para prever os impactos causados, por exemplo, pela difusão
de um determinado contaminante, em um aquífero ou pela explotação desordenada
da água subterrânea. Para tanto, são vários os programas como o MATLAB,
SCILAB, HEC HMS, MOD FLOW, que auxiliam na modelagem do problema.
2.3.1 Propagação de poluentes no meio aquático
Os poluentes ao atingirem o meio aquoso não se comportam, em geral,
de maneira estática sendo, portanto sujeitos a uma série de fatores, físicos,
químicos, e biológicos existentes na natureza, os quais alteram sua composição e
concentração.
Entre os fenômenos de transporte, podemos destacar a própria diluição
que é uma propriedade característica dos líquidos que em via de regra reduz a
concentração original do poluente e os efeitos hidrodinâmicos.
7
A ação hidrodinâmica, é realizada pela característica intrínseca aos
líquidos de transporte de poluente de seu despejo a outras regiões e, por
conseguinte ocasionando uma variação na sua concentração no espaço e no tempo.
Sendo o carreamento desses poluentes feitos pelo campo de velocidades da água,
ou como é mais conhecido pelo processo da advecção. Além disso, processos de
concentração de substâncias; quer sejam dissolvidas ou em suspensão é devido a
um fenômeno denominado de difusão. Estes, por sua vez, de acordo com Braga
(2005) dividem-se em dois tipos básicos saber:
1. Difusão Molecular: É um fenômeno independente de forças externas e
que consiste no movimento característico das moléculas do fluido, por
conta da agitação térmica de suas partículas e que geralmente em
aquíferos naturais não causa grande impacto em sua concentração.
2. Difusão Turbulenta: É característico de um escoamento turbulento da
água, ou seja, com a intensificação no gradiente de velocidade do líquido
pode ocorrer um espalhamento maior no poluente, uma diminuição em sua
concentração e promoção de uma aeração superior e, por conseguinte um
aumento no Oxigênio Dissolvido (OD) e queda da Demanda Bioquímica de
Oxigênio (DBO). Além de uma diminuição no efeito nos gases produzidos
por organismos em decomposição anaeróbica, como podemos citar o gás
metano e o ácido sulfídrico.
Contudo, além dos processos supracitados, existe um fenômeno de
propagação e transporte de poluentes causados pela ocorrência conjunta de difusão
molecular, concomitante com a turbulenta e da própria advecção, denominado de
dispersão.
Segundo Braga (2005) um dos primeiros modelos matemáticos sobre
qualidade de água é denominado de Modelo de Streeter-Phelps. Este modelo é
consistente com o fenômeno de decaimento radioativo, meia vida de uma
substância, crescimento populacional, ou seja, todos esses fenômenos têm uma
taxa de variação, para mais ou para menos, proporcional à quantidade existente em
um dado instante, isto é, a matéria orgânica biodegradável tem uma taxa de
decomposição proporcional à concentração de matéria orgânica em um dado
8
instante de tempo. Então assim, a descrição do modelo representa uma equação
linear de primeira ordem separável.
(2.1)
Onde:
→ é a demanda bioquímica de oxigênio;
→ é a constante de desoxigenação à qual depende do efluente;
Separação de variáveis:
(2.2)
Integra-se em t
∫
∫
(2.3)
Gera uma função logarítmica
(2.4)
A DBO tem uma taxa de decaimento com modelo do tipo função exponencial
(2.5)
9
2.4 Equações diferenciais
Segundo Zill (2001, p. 2) a definição de equação diferencial, consiste em
uma equação que contêm derivadas, ou diferenciais de uma ou mais variáveis
dependentes, em relação a uma ou mais variáveis independentes. As equações
diferenciais são classificadas pela ordem, tipo e linearidade.
2.4.1 Classificação de equações diferenciais quanto ao tipo
Se uma equação contém somente derivadas ordinárias, com relação a
uma única variável dependente, a mesma é dita Equação Diferencial Ordinária
(EDO), como exemplo:
(2.6)
(2.7)
Uma equação que envolve as derivadas parciais de uma ou mais variáveis
dependentes de duas ou mais variáveis independentes é chamada de equação
diferencial parcial (EDP), como exemplo:
(2.8)
(2.9)
2.4.2 Classificação de equações diferenciais quanto à ordem
A ordem da equação diferencial é a ordem da derivada de mais alta
ordem da equação. Sendo assim, as equações (2.6), (2.7), (2.8) e (2.9), são de
ordem 2, 4, 2 e 4, respectivamente.
10
2.4.3 Classificação de equações diferenciais quanto à linearidade
Uma equação diferencial ordinária é chamada de linear quando pode ser
escrita na forma:
(2.10)
Observe que as equações diferenciais lineares são caracterizadas por
duas propriedades:
(i) A variável dependente e todas as suas derivadas são do primeiro grau;
(ii) Cada coeficiente depende apenas da variável independente
2.4.4 Solução de equação diferencial linear de 1ª ordem
Uma equação linear de primeira ordem é toda equação do tipo:
(2.11)
É importante sabermos solucionar essa equação, haja vista muito
problemas serem modelados por ela, como exemplo podemos citar, problemas de:
decaimento radioativo, crescimento de bactérias, mistura de líquidos, intrusão salina,
datação de carbono 14, entre outros. Então assim, para determinamos uma
solução para (2.11), vamos utilizar uma função auxiliar , por enquanto
desconhecida, à qual é comumente chamada de fator integrante ou fator de
integração da equação linear. Multiplicando, então a equação (2.11) por e
impondo que o resultado à esquerda da igualdade, seja a derivada com relação a
do produto de . Temos que:
(2.12)
Então,
(2.13)
11
Assumindo e integrando em ,
∫
∫ ∫
(2.14)
Tomando e usando a definição de logaritmo, determinamos a função auxiliar,
∫
(2.15)
Supondo (2.15) voltamos para a equação (2.12) e conclui-se que:
(2.16)
Integrando, chega-se que:
∫
(2.17)
E, enfim, chegamos à solução geral:
∫
∫
(2.18)
2.4.5 Equações diferenciais de 2ª ordem
Para solucionarmos uma equação diferencial de 2ª ordem, do tipo linear e
não homogênea, faremos uso do Método da Variação de Parâmetros (MVP).
Seja a equação de segunda ordem
+ + (2.19)
Admitindo , colocamos na forma dita padrão:
+ + (2.20)
Supondo que e sejam solução de (2.20), sabe-se
que a combinação linear, abaixo, em que e , são constantes denominadas
parâmetros, também é solução de (2.20), à qual recebe o nome de função
complementar:
(2.21)
12
Chamaremos de , uma solução particular de (2.20), que obteremos
variando os parâmetros e , isto é:
(2.22)
Para substituir (2.22) em (2.20), precisamos obter e
, daí, impondo
para simplificar os cálculos de primeira e segunda derivada:
(2.23)
Derivando (2.22) e usando (2.23), temos:
+
(2.24)
Com isso,
+
(2.25)
Substituindo (2.25) em (2.20), segue-se que:
+
+
[
] +
+
(2.26)
Mas como e são soluções de (2.20), tem-se:
+
+ +
+ (2.27) Dessa forma, chegamos ao seguinte sistema de Cramer, nas variáveis
e :
,
(2.28)
A solução em questão pode ser expressa através de determinantes:
;
(2.29)
Onde
|
|, |
| e |
|, (2.30)
13
O determinante é chamado de Wronskiano e quando é não nulo indica
que o conjunto { é Linearmente Independente (LI) e, portanto, é um conjunto
solução da EDO. Com isso, a solução geral de (2.20), é a soma da função
complementar , com a solução particular .
(2.31)
2.4.6 Equações diferenciais de ordem superior
Podemos generalizar o método anterior para equações lineares
homogêneas de ordem n, da seguinte forma:
(2.32)
Como a ordem é de ordem n a função complementar é a combinação
linear, abaixo:
(2.33)
Com isso, temos uma solução particular para (2.32) é
(2.34)
Podemos determinar os , por uma sistema de Cramer de
n equações e n incógnitas. De modo que o cálculo dos , é feito naturalmente
através de integral na variável , então temos que:
{
(2.35)
Com isso, pela regra de Cramer, concluímos que:
(2.36)
14
Sendo que é o Wronskiano e é o Wronskiano, cuja i-ésima
coluna tem como coeficientes:
(2.37)
2.4.7 Aproximação em série de potências para uma equação diferencial
linear de segunda ordem
Para se encontrar a solução geral de uma equação diferencial de 2ª
ordem faz-se necessário à determinação de um conjunto fundamental de solução, à
qual geralmente é formada por funções elementares analíticas linearmente
independentes, ou seja, o Wronskiano do conjunto fundamental deve ser não nulo e,
cada função, pode ser expressa em uma série de potências, ou série de Taylor na
vizinhança de um determinado ponto.
Em problemas, envolvendo Equações Diferenciais de 2ªordem com
coeficientes constantes, podemos utilizar o método dos coeficientes a determinar ou,
ainda, o método das variações dos parâmetros. Contudo, em muitos casos práticos
e problemas de engenharia, tais como de estrutura, deflexão de vigas,
amortecimento e de hidráulica, as equações não costumam ter equações com
coeficientes constantes e, por conseguinte não comumente apresentam como
soluções as funções elementares conhecidas. Dessa forma, é importante o
desenvolvimento em séries para aproximar uma solução analítica, para tanto é vital
que saibamos a convergência ou divergência de uma determinada série.
15
2.4.8 Série de potências
Definimos série de potências em é uma série da forma:
∑
(2.38)
Uma série de potências converge em um determinado ponto , se o limite
da série existir e for finito. Isto é:
∑
(2.39)
É possível mostrar que toda série absolutamente convergente é
convergente, mas a recíproca não é verdadeira.
∑| |
∑| || |
(2.40)
Então,
∑
(2.41)
Teste da razão para a convergência em módulo:
Se e se para um valor fixo de :
|
| | | |
| | | (2.42)
Então a série em questão é convergente, caso | | ; divergente se
| | e inconclusivo quando | | .
16
2.4.8.1 Raio de convergência de uma série de potência
Existe um número não negativo , chamado de raio de convergência, tal
que ∑
, converge absolutamente para | | e diverge para
| | . Caso a série convirja em apenas um ponto , define-se ;
caso a série convirja para todo , dizemos que é infinito. Então assim, se
tivermos , o intervalo | | é conhecido como intervalo de convergência,
ou vizinhança do ponto .
Figura 2.2 – Raio de convergência de uma série
Fonte: Figura modificada, Boyce (2010, p. 132.)
Encontrar uma solução para EDO em série de potência:
(2.43)
A solução é o conjunto independente , ou seja a combinação
linear dos elementos desse conjunto linearmente independente, cujo wronskiano é
diferente de zero, a saber:
(2.44)
17
Esse tipo de equação não esboça dificuldades, mas serve para mostrar a
aplicação do método em série. Isto é, tome
∑
(2.45)
Derivando termo a termo temos que:
∑
(2.46)
Derivando novamente,
∑
(2.47)
Substituindo na equação original
∑
∑
(2.48)
Trocamos o índice do primeiro somatório substituindo por ,
iniciando agora do zero ao invés de dois, segue-se que:
∑
(2.49)
Para que a soma seja nula é necessário e suficiente que
(2.50)
Temos então uma relação de recorrência para índices pares, na qual variamos os
valores de :
(2.51)
(2.52)
18
(2.53)
(2.54)
Para os índices ímpares, temos uma relação análoga:
(2.55)
(2.56)
(2.57)
(2.58)
A rigor teríamos que provar a validade das fórmulas de recorrência por
meio de indução matemática. Contudo, por se tratar de um requisito básico e
elementar, aceitaremos sem demonstrar esse fato.
Substituindo os coeficientes encontrados na fórmula de série de potências:
∑
(2.59)
Por conseguinte:
(2.60)
+...
19
O desenvolvimento de uma função em série de Taylor ou Maclaurin
na vizinhança do ponto , é dado por:
∑
(2.62)
É fácil ver que a expansão em série de Taylor da função é:
(2.63)
A expansão em série de Taylor da função :
(2.64)
Conclui-se que a solução da equação diferencial de 2ª ordem:
(2.65)
É dada pela combinação linear das funções e de , assim:
(2.66)
*
++ *
+
(2.61)
20
2.4.9 Série de Fourier
Muitos problemas que envolvem Equações Diferenciais Parciais (EDP) são
solucionados com o uso de séries infinitas. Entre elas a conhecida série de Taylor, e
as que são combinações de senos e cossenos e são chamadas de séries de
Fourier, às quais, exercem um papel importante no método de separação de
variáveis e são escritas nessa forma:
∑ (
)
(2.67)
As séries de Fourier também são vistas em uma ampla gama de
problemas físicos e de Engenharias, entre os quais, temos os de análises de
sistemas mecânicos e de molas, bem como sistemas elétricos.
Uma propriedade bem conhecida de funções seno e cosseno é da sua
periodicidade de radiano, isto é de maneira geral uma função é dita periódica
com período se o domínio de contém , sempre que contiver e se:
(2.68)
Outra propriedade importante de funções seno e cosseno é a sua
ortogonalidade dessa forma escreve-se um produto dito interno ou escalar padrão
de duas funções reais no intervalo , como:
∫
(2.69)
Observa-se que analogamente ao caso de vetores o produto escalar é dito
ortogonal se , isto é:
∫
(2.70)
As funções
e
; ; seguem as seguintes relações
de ortogonalidade:
I. Caso ,
∫
(2.71)
21
II. Caso ,
∫
(2.72)
Para quaisquer e ,
∫
(2.73)
III. Caso ,
∫
(2.74)
IV. Caso ,
∫
(2.75)
2.4.10 Equações diferenciais parciais e o método de separação de variáveis.
Em muitos ramos da Física, da Engenharia e da própria Matemática, nos
quais são envolvidos fenômenos como difusão, propagação do calor, processos
oscilatórios e teoria do potencial. Assim como, outros problemas estacionários e
transientes envolvem, diretamente, equações diferenciais parciais, entre as quais, a
equação do calor, a equação de Laplace e a equação da onda.
Um dos métodos mais antigos para solucionar equações desse tipo, é o
Método da Separação de Variáveis. Este método foi usado por volta do ano de 1750,
por vários matemáticos de renome como D’Alembert, Daniel Bernoulli e Euler
(BOYCE, 2002, p. 311).
Para solucionarmos o clássico problema da propagação de calor por uma
barra sólida, utilizaremos o método da separação de variáveis e recairemos em uma
série de Fourier. Série esta, é semelhante a serie de Taylor no sentindo de também
22
ser usada para cálculos de equações diferenciais, através de aproximação
numérica. Contudo, ao invés de ser uma série que combina potências com as
derivadas da função no ponto de interesse, a mesma aproxima harmonicamente,
através da combinação de senos e cossenos.
Seja, abaixo, a Figura 2.3, que ilustra uma barra, na qual há uma
propagação de calor, evidenciando a função , dependente da distância e do tempo
e suas respectivas, condições de contorno estacionárias nas extremidades.
Figura 2.3 – Barra de propagação de calor
Fonte: Figura Modificada, Boyce (2002, p. 311)
A equação do calor tem a forma:
(2.76)
Onde é uma constante conhecida como difusividade térmica, à qual
depende da natureza do material do qual a barra é feita e é definido por:
(2.77)
Sendo , denominado condutividade térmica, é a densidade e é o calor
específico do material na barra.
23
Para resolvermos o problema temos que encontrar uma função do tipo
que satisfaça (2.76) e que obedeça a condição inicial, abaixo, com uma
função dada por:
(2.78)
Suponha ainda que a temperatura seja constante igual à zero nas
extremidades, é o que chamamos de condições de contorno, a saber:
(2.79)
Para usarmos o método conhecido como separação de variáveis, temos
que primeiro supor que pode ser escrita como produto:
(2.80)
Assim substituímos (2.80) em (2.76) e obtemos:
(2.81)
Separando as variáveis de , tem-se
(2.82)
Igualamos a uma constante – , chamada de constante de separação,
(2.83)
Reduzimos o problema de EDP para EDO:
{
(2.84)
Uma solução para a equação (I) do sistema (2.84) é encontrada utilizando
uma equação característica:
(2.85)
24
Ou seja,
√ √ (2.86)
Sendo assim, é sabido que:
√ √ (2.87)
Em face disso, usando a condição de contorno,
(2.88) √
(2.89)
Como , e √ , logo como a função seno se anula para os
múltiplos inteiros de . Conclui-se que:
√
(2.90)
Substituindo em (2.84) item (II), temos,
(2.91)
Resolvendo a equação, chega-se que:
(2.92)
Usando o fato expresso em (2.80) chega-se que,
(2.93)
Contudo, ainda resta obedecer à condição (2.78), isto é, temos que formar uma série
da combinação linear do conjunto solução da equação de tal forma que os
coeficientes satisfaçam às condições iniciais do problema em questão.
25
Vamos tomar então,
∑
(2.94)
E igualando a
∑
∑ (
) (
)
(2.95)
Evidentemente temos que admitir que a série convirja e verifique a
condição inicial do problema posto, com isso,
∑ (
)
(2.96)
A série encontrada é a serie de Fourier, que está convergindo para a
condição inicial de distribuição de temperatura . Cujos coeficientes são
devidamente encontrados e estão bem definidos na forma:
∫
(2.97)
Tabela 2.4 – Valores de Difusividade Térmica
Fonte: Boyce (2002, p. 311)
VALORES DE DIFUSIVIDADE TÉRMICA PARA ALGUNS MATERIAIS COMUNS
Material (
Prata 1,71
Cobre 1,14
Alumínio 0,86
Ferro fundido 0,12
Granito 0,011
Tijolo 0,0038
Água 0,00144
26
3. METODOLOGIA
Os métodos numéricos são muito importantes, haja vista à maioria dos
problemas que surgem em Física e em Engenharia consistem de Equações
Diferenciais. Quer sejam Equações Diferenciais Ordinárias ou Equações Diferenciais
Parciais, muitas têm um nível de complexidade bem alto de tal modo que é
extremamente difícil, ou em muitos casos impossível, encontrarmos uma solução
analítica para o problema.
Existem vários métodos numéricos que servem de ferramenta para
resolver as referidas equações diferenciais, entre eles, existem: diferenças finitas,
elementos finitos, elementos de contorno e, agora mais recentemente, tem recebido
a devida atenção de estudiosos, o método conhecido como “sem malhas”
(MESHLESS). Isto é devido ao fato de não ser necessário a discretização dos
pontos do domínio da função, de modo a formar malhas e também, por conta de
suas propriedades de convergência.
O método que iremos utilizar nesse trabalho é o método de Kansa, no qual
consiste em uma abordagem, relativamente, nova, por volta de 1990, para solução
de equações diferenciais parciais através das funções de base radial, RBF (Radial
Basis Functions). Nesse ínterim, modelaremos desde uma equação diferencial
ordinária linear de primeira ordem, até a equação diferencial parcial transiente, que
descreve a concentração de um poluente conservativo, em função do tempo e do
espaço percorrido, denominada de equação da difusão-advecção.
27
3.1 O método sem malhas: Funções de Base Radial
Uma função RBF é chamada assim por conta da sua dependência, ou
simetria radial, da distância entre um ponto central da função e um ponto genérico
denotado por .
As funções de base radial se apresentam da seguinte forma:
( ) (‖ ‖) (3.1)
A norma ‖ ‖ representa a distância entre e , à qual chamaremos
de r. As funções de base radial a serem estudadas também dependem de um
parâmetro de forma, c, sendo então expressas da seguinte maneira:
(‖ ‖ ) (3.2)
Onde, respectivamente o caso unidimensional e o n-dimensional:
[
]
(3.3)
*
( )
+
(3.4)
Alguns exemplos de Funções de Base Radial (RBF) estão dispostos na tabela a
seguir:
Tabela 3.1 – Funções de base radial
RBF Tipo de função
Splines de Placas Finas | | | | , com par
Multiquadráticas √
Multiquadráticas Inversas
Gaussianas
Fonte: Liu (2003)
28
De acordo com Ferreira (2009) o método trabalha com pontos espalhados
no domínio de interesse e o interpolador RBF é uma combinação linear de funções
de base radial centrada nos pontos :
∑
(3.5)
Sendo que é número de pontos de dados; os
são os coeficientes a serem determinados; e é a função de base radial.
Um problema a ser transposto é determinar um , parâmetro de forma
ótimo, de modo a reduzir ao máximo o erro e, por conseguinte o resíduo na
discretização de pontos no computador, haja vista linguagem de computação
entender apenas a Aritmética discreta. Contudo, trataremos desse problema de
determinação do parâmetro de forma, a posteriori.
3.2 Operadores lineares aplicados em uma EDP
Seja um problema do tipo EDP elíptica:
com (3.6)
Com condições de contorno do tipo:
com (3.7)
Onde:
e são operadores lineares que dependem do tipo de problema a ser estudado.
é um operador linear parcial diferencial. É importante lembrar que
.
Por exemplo, para uma equação do tipo:
(3.8)
29
Pode-se escrever esta equação na forma, , onde:
(3.9)
Se a condição de contorno for do tipo Dirichlet então,
(3.10)
3.3 Método de Kansa
Segundo Ferreira (2009) as primeiras tentativas de resolver EDP com
técnicas RBF foram feitas por Kansa, por volta de 1990. Desde então, vários
estudos foram desenvolvidos para investigar a aplicabilidade deste método para
solução numérica de EDP, no qual, de acordo com Liu (2001) e Ferreira (2009),
utilizou-se a seguinte função de aproximação:
∑
(3.11)
Para um problema do tipo EDP elíptica, estacionário, o qual independe do tempo:
; (3.12)
Com condições de contorno de um problema de Dirichlet:
; (3.13)
onde e são operadores lineares e é a fronteira, ou contorno de .
Substituindo a função de aproximação na equação diferencial, tem-se:
[∑
]
(3.14)
Sendo L um funcional linear,
[∑
]
(3.15)
30
Substituindo a função de aproximação na condição de contorno
∑
(3.16)
Ou, simplesmente:
∑
(3.17)
Assim, o problema se resume a resolver um sistema na forma matricial do tipo:
(3.18)
A matriz de colocação A é formada por dois blocos:
[
] (3.19)
Os dois blocos são gerados da seguinte maneira:
[ (‖ ‖ )] tal que (3.20)
(‖ ‖ ), com, (3.21)
O conjunto X é dividido em pontos do contorno e pontos do interior de Ω. O
vetor y é constituído das entradas de e , da seguinte forma:
[
] (3.22)
Pode-se notar que uma mudança nas condições de contorno da matriz
acarreta em uma simples mudança de algumas linhas da matriz de colocação A e do
vetor do lado direito, y. De forma mais clara, tem-se um sistema matricial da seguinte
forma:
[
] [ ] [
] (3.23)
31
O sistema é então resolvido para os e esses valores são em
seguida substituídos na função de aproximação.
(3.24)
Qualquer equação diferencial, com ou sem solução analítica conhecida
pode ser, solucionada com uso do Método numérico sem malhas. Nessa
perspectiva, a posteriori, trataremos primeiro de uma equação diferencial ordinária
linear, de 1ª ordem, depois evoluindo para de 2ª ordem e por fim, faremos uso de um
problema usando EDP, de 2ª ordem transiente. Para tanto, faz-se necessário o uso
de uma ferramenta computacional, como o SCILAB, para desenvolvermos uma
rotina de otimização do parâmetro de forma c.
3.4 Equação geral apresentada no SCILAB
(
) (
) (3.25)
Inicialmente implementamos a rotina que solicita ao usuário os seguintes
parâmetros de entrada:
Números de pontos do domínio;
Os limites inferior e superior do domínio;
Valores das condições de contorno;
Parâmetro de forma “ ”.
Atribuir em seguida os valores das constantes envolvidas na EDO, ou seja,
os valores de e . De maneira a deixar claro para o programa a
forma da EDO envolvida no problema.
Os passos seguintes são as funções internas do programa para calcular
na sequência:
32
a) Os vetores com os pontos do domínio;
b) Determinar a EDO envolvida;
c) Determinar a matriz de distâncias ;
d) Determinar a matriz :
√ (3.26)
Onde,
√( )
(3.27)
e) Obter a matriz ), que varia em função da EDO envolvida:
(3.28)
f) Montar a matriz A, a partir da matriz e :
nn,n,1
n1,n1,1n
n2,2,1
n1,1,1
Φ...Φ
:...:
ΦL...ΦL
ΦL...ΦL
Φ...Φ
(3.29)
g) Em seguida o programa verifica se a matriz formada possui inversa e a
gera, se sim ele prossegue, caso contrário envia uma mensagem informando
que a matriz gerada não possui inversa;
h) Determina o vetor y;
i) Determina os valores de λ;
j) Gera a matriz de aproximação.
(3.30)
k) Plota os resultados em gráfico.
33
Seguem, doravante, as legendas com as quais denotamos os elementos e
matrizes do sistema:
e → São as matrizes que armazenam os pontos envolvidos na discretização
dos pontos do domínio e j;
→ Valor da distância entre os pontos e ;
→ Matriz com os valores das distâncias ;
→ Contêm os valores da matriz ;
→ Contêm os valores da matriz L(Ф);
→ Contêm os valores da matriz L(Ф) para o termo da 1ª derivada;
→ Contêm os valores da matriz para o termo da 2ª derivada;
→ Contêm os valores da matriz A;
→ Contêm os valores da matriz inversa de A;
→ Contêm os valores do vetor y;
→ Contêm os valores dos multiplicadores encontrados;
→ Contêm os valores da matriz de aproximação;
34
4. SIMULAÇÕES E RESULTADOS
Utilizaremos primeiro, para ilustrar, o método sem malhas (MESHLESS)
em uma Equação Diferencial Linear de 1ª Ordem, EDO. De início a solução da EDO,
abaixo, é simples de ser obtida utilizando o método do fator de integração. Em
seguida, utiliza-se o método sem malhas, propriamente dito, com a função de base
radial, mais conhecida como método de Kansa. Então gera-se as matrizes do
problema, plota-se em um gráfico a solução numérica e a solução analítica com o
objetivo de observarmos o comportamento e o ajustamento das curvas, concluindo
portanto que nesse caso a suposição do parâmetro de forma c, como sendo a razão
entre dois e a raiz quadrada de n, onde este é o número de pontos, serviu nesse
primeiro momento para estimarmos um valor próximo a solução real do problema
apresentado.
Doravante, utilizaremos um problema mais complexo utilizando a Equação
de Difusão Pura e da Equação de Difusão Advecção, à qual separaremos as
variáveis e encontraremos a sua solução quase analítica, via série de Fourier, para
por fim usarmos de diferenças finitas e aplicarmos a rotina MESHLESS, ou método,
sem malhas, no SCILAB, para modelarmos o problema, apresentarmos as suas
matrizes e plotarmos os resultados da comparação entre a solução analítica e a
solução numérica, via método de Kansa.
4.1 Aplicação do método sem malhas (MESHLESS) em uma EDO via Método de
Kansa.
Nesse exemplo, vamos aplicar o método sem malhas com função de base
radial (RBF) e com a aproximação de Kansa, para resolver uma equação diferencial
ordinária e de primeira ordem.
Iniciaremos a aplicação do método numérico sem malhas, aplicando-o a
uma Equação Diferencial Linear de 1ª ordem, EDO, descrita a seguir, com sua
respectiva condição de contorno:
(4.1)
35
É fácil mostrar que a solução da equação diferencial é:
(4.2)
Nesse caso, em particular, a solução analítica é conhecida, mas com o
desenvolvimento de estudos através dos métodos sem malha, é possível estimar
uma solução numérica, mesmo sem conhecer a solução analítica. Neste trabalho,
porém, nossa missão é de comparar a solução numérica, obtida através de métodos
sem malha, com a solução analítica de algumas equações diferenciais. O estudo
requer a otimização do parâmetro de forma, de modo a melhorar a precisão e, por
conseguinte, tornar o erro pequeno. Nessa perspectiva, o princípio dos métodos sem
malha baseia-se da não necessidade de utilizarmos uma relação entre os nós, assim
como o próprio nome sugere, evitaremos o traçado de uma malha na discretização
dos pontos. Mas o que vem a ser uma malha? Segundo LIU: “Uma malha é definida
como qualquer um dos espaços abertos ou interstícios entre os fios de uma rede
que é formada ligando os nós de uma forma pré-definida” (LIU, 2003, p.19).
Para a aplicação do método sem malhas com aproximação de Kansa,
considere a seguinte solução aproximada ( ):
∑ ( )
(4.3)
Sendo que:
O número de pontos ( ;
Linha ;
Coluna ;
Parâmetro de forma ( ;
Matriz das distâncias com N linhas e N colunas
Matriz da função de base radial com N linhas e N colunas (
Vetor com N elementos que contém os coeficientes da solução aproximada
( ) √
(4.4)
36
( ) (4.5)
Para cada ponto ( ):
{
( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( )
(4.6)
N
1
NN,N,1
N1,-N1,1-N
N2,2,1
N1,1,1
1
2
1
c),(rc),(r
c),(rc),(r
c),(rc),(r
c),(rc),(r
)(
)(
)(
)(
N
N
xY
xY
xY
xY
(4.7)
Na forma compacta:
(4.8)
Para aplicar a aproximação de Kansa ao problema, observe os passos
para a resolução:
1. Escolher um intervalo para a resolução da EDO. Nesse caso, vamos trabalhar com
um intervalo [0,2].
2. Escolher o número de pontos e a sua localização. Considere inicialmente 7 pontos,
distribuídos em intervalos iguais.
3. O valor do parâmetro de forma é arbitrário, nesse caso, adotaremos a seguinte
aproximação (FERREIRA, 2009 p.87):
√
√
(4.9)
No qual, N é o número de pontos.
37
4. Aplicar a aproximação de Kansa
Observe que o vetor Lambda, que representa o coeficiente para a aproximação,
apresenta um valor constante, considerando um parâmetro de forma conhecido,
Portanto:
(4.10)
Substituindo na equação diferencial:
(4.11)
Organizando a equação:
*
+
(4.12)
A equação pode ser reescrita da seguinte forma:
(4.13)
Tem-se que:
(4.14)
(4.15)
(4.16)
5. Desenvolver os termos da EDO
A) Desenvolvendo para FI:
( ) (4.17)
38
Desenvolvendo a equação:
√( ) √
(4.18)
B) Desenvolvendo para DFI:
Observamos que é a variável , assim a derivada de FI pode ser expressa como:
[ ( )]
√( )
(4.19)
Organizando a equação:
[( )]
√( )
√
(4.20)
6. Substituir na equação diferencial:
√
√
(4.21)
Substituindo na equação diferencial:
*√
√
+ (4.22)
Na forma matricial, temos que:
N
1
NN,N,1
N1,-N1,1-N
N2,2,1
N1,1,1
LFILFI
LFILFI
LFILFI
LFILFI
)(
))1(
)2(
)1(
NX
NX
X
X
(4.23)
39
7. Condição de contorno
Nesse exemplo, e EDO apresenta a condição de contorno . Com base na
aproximação de Kansa, a solução aproximada pode ser expressa por:
( ) ( ) ( ) (4.24)
Na forma compacta:
(4.25)
8. Agrupando as condições de contorno com a EDO, temos que para a determinação
dos coeficientes do vetor LÂMBDA, faz-se preciso a resolução de um sistema de
equações. Nesse sentindo, agruparemos as equações, substituindo a primeira linha
da matriz da EDO pela sua condição de contorno. Feito isso, segue-se que:
{
( ) ( ) ( )
(4.26)
Ou, equivalentemente, em sua forma matricial, temos que:
N
1
NN,N,1
N1,-N1,1-N
N2,2,1
N1,1,1
LFILFI
LFILFI
LFILFI
FIFI
)(
))1(
)2(
)1(
NX
NX
X
Y
(4.27)
Observe que, nesse caso, a condição de contorno com e
corresponde a primeira linha da matriz. Na forma compacta:
(4.28)
40
Nesse caso:
A= Matriz com os coeficientes FI (primeira linha) e LFI (demais linhas)
F= Vetor com o coeficiente Y(1), para a primeira linha e X(i) para as demais linhas,
tem-se X(i) representa o vetor para a variável x, no intervalo [0,2].
9. Montar as matrizes e Vetores (para N=7)
Vetores , , dispostos na tabela 4.1, abaixo:
Tabela 4.1 – Discretização dos pontos do domínio da EDO de 1ª ordem
1 2 3 4 5 6 7
1 0 0.3333 0.6667 1.0000 1.3333 1.6667 2.0000
2 0.3333
3 0.6667
4 1.0000
5 1.3333
6 1.6667
7 2.0000
Matriz R= :
( ) (4.29)
Tabela 4.1 – Matriz R das distâncias
1 2 3 4 5 6 7
0,000000 -0,333333 -0,666667 -1,000000 -1,333333 -1,666667 -2,000000
0,333333 0,000000 -0,333333 -0,666667 -1,000000 -1,333333 -1,666667
0,666667 0,333333 0,000000 -0,333333 -0,666667 -1,000000 -1,333333
1,000000 0,666667 0,333333 0,000000 -0,333333 -0,666667 -1,000000
1,333333 1,000000 0,666667 0,333333 0,000000 -0,333333 -0,666667
1,666667 1,333333 1,000000 0,666667 0,333333 0,000000 -0,333333
2,000000 1,666667 1,333333 1,000000 0,666667 0,333333 0,000000
41
Matriz FI = :
√
(4.30)
Tabela 4.11 – Matriz FI = ϕ, função de interpolação, ponto a ponto.
1 2 3 4 5 6 7
0.755929 0.826146 1.007927 1.253566 1.532683 1.830114 2.13809
0.826146 0.755929 0.826186 1.007927 1.253566 1.53277 1.830114
1.007927 0.826186 0.755929 0.826146 1.007861 1.253566 1.532683
1.253566 1.007927 0.826146 0.755929 0.826146 1.007927 1.253566
1.532683 1.253566 1.007861 0.826146 0.755929 0.826186 1.007927
1.830114 1.53277 1.253566 1.007927 0.826186 0.755929 0.826146
2.13809 1.830114 1.532683 1.253566 1.007927 0.826146 0.755929
Matriz DFI:
√
(4.31)
Tabela 4.12– Matriz DFI das derivadas da função FI, ponto a ponto.
1
0 -0,40347 -0,66144 -0,79772 -0,86992 -0,91071 -0,93541
0,403473 -0,40347 -0,66144 -0,79772 -0,86992 -0,91071
0,661438 0,403473 0 -0,40347 -0,66144 -0,79772 -0,86992
0,797724 0,661438 0,403473 0 -0,40347 -0,66144 -0,79772
0,869918 0,797724 0,661438 0,403473 0 -0,40347 -0,66144
0,910705 0,869918 0,797724 0,661438 0,403473 0 -0,40347
0,935414 0,910705 0,869918 0,797724 0,661438 0,403473 0
42
Matriz LFI
(4.32)
1 2 3 4 5 6 7
0.755929 0.422707 0.346471 0.455842 0.662771 0.919406 1.202676
1.229586 0.755929 0.422646 0.346471 0.455842 0.662842 0.919406
1.669384 1.229727 0.755929 0.422707 0.346461 0.455842 0.662771
2.05129 1.669384 1.229586 0.755929 0.422707 0.346471 0.455842
2.402595 2.05129 1.669262 1.229586 0.755929 0.422646 0.346471
2.402698 2.05129 1.669384 1.229727 0.755929 0.422707
2.740822 2.402595 2.05129 1.669384 1.229586 0.755929
Matriz A:
(4.33)
1 2 3 4 5 6 7
0.755929 0.826146 1.007927 1.253566 1.532683 1.830114 2.13809
1.229586 0.755929 0.422646 0.346471 0.455842 0.662842 0.919406
1.669384 1.229727 0.755929 0.422707 0.346461 0.455842 0.662771
2.05129 1.669384 1.229586 0.755929 0.422707 0.346471 0.455842
2.402595 2.05129 1.669262 1.229586 0.755929 0.422646 0.346471
2.740822 2.402698 2.05129 1.669384 1.229727 0.755929 0.422707
3.073504 2.740822 2.402595 2.05129 1.669384 1.229586 0.755929
Tabela 4.13 – Matriz LFI, soma das matrizes FI e DFI
Tabela 4.14– Matriz de colocação, A, é o produto da inversa de λ pela matriz F
43
10. Inverter a matriz A, determinaremos o coeficiente LAMBDA e a solução numérica.
Matriz :
(4.34)
1 2 3 4 5 6 7
-0.73976 8.641945 -15.5759 17.63564 -16.5333 13.48034 -5.35704
0.072258 -18.016 40.85807 -52.1484 50.94688 -42.11 17.52807
1.947293 12.09544 -40.0311 66.23726 -72.7824 62.59557 -26.7075
-2.65666 -3.34656 19.36612 -49.0664 69.92078 -67.9038 30.11678
2.974397 -0.1808 -4.61464 21.72549 -45.8846 58.43369 -28.8928
-2.61038 1.201901 -0.63025 -5.01575 16.9383 -33.4516 20.4409
1.443143 -0.73317 1.083888 0.176638 -2.13043 8.553822 -6.73079
-3.1174 0.0000 4.0000 4.0000 0
2.1742 0.3333 2.9160 3.0171 0,101087
2411 0.6667 2.2338 2.2934 0,05964
-7.6124 1.0000 1.8394 1.8898 0,050411
8.9147 1.3333 1.6513 1.6805 0,029199
-7.7657 1.6667 1.6110 1.6395 0,028501
4.3817 2.0000 1.6767 1.6860 0,009293
Tabela 4.15– Inversa da matriz de colocação A
Tabela 4.16 – Resumo dos dados da solução da EDO de 1ª ordem via método sem malhas
44
4.1.1 Vetor Resíduo - VR
O vetor resíduo é definido pelo produto da matriz A pela matriz ,
subtraído pela matriz F. Sua importância deve-se ao fato de servir como indicador da
qualidade da aproximação entre a solução numérica e a solução analítica de uma
determinada equação diferencial. Isto por que, para cada c, parâmetro de forma,
gera-se um vetor resíduo e quanto menor for o somatório do módulo das
componentes deste, melhor será a aproximação e, então, podemos identificar o
parâmetro de forma c – ótimo. Nesse exemplo, em particular, o parâmetro de forma
foi arbitrário e igual a 0,755929 (vide item 4.9), contudo esses valores sugeridos na
literatura nem sempre refletem em um resultado razoável, por isso, a posteriori,
faremos o uso da plataforma SCILAB 5.4.1, para otimizarmos o parâmetro de forma.
Então, nesse caso particular, a somatória dos módulos das componentes
do vetor resíduo , é o número ínfimo, apresentado abaixo, que caracteriza uma
aproximação significativamente razoável:
∑| |
(4.35)
Eis abaixo, nas figuras, 4.1; 4.2 e 4.3, a comparação gráfica entre a
solução numérica e a solução analítica, através da simulação com, respectivamente,
5, 7 e 10 pontos de discretização do domínio, no qual a equação está definida e com
um parâmetro de forma c suposto igual a 0,755929.
Tabela 4.17– Componentes do vetor resíduo, VR, da EDO de 1ª ordem
= (A
-7,10543x10-15
-9,4369x10-16
2,18714x10-14
1,64313x10-14
6,88338x10-15
5,77316x10-15
9,32587x10-15
45
Figura 4.1 – Parâmetro de forma utilizado "c=0,755929" , pontos do domínio, com número de pontos N=5.
46
Figura 4.2 – Parâmetro de forma utilizado "c=0,755929" , pontos do domínio, com número de pontos N=7.
Figura 4.3– Parâmetro de forma utilizado "c=0,755929" , pontos do domínio, com número de pontos N=10.
47
4.2 Aplicação do método sem malhas (MESHLESS) em uma EDO de 2ª ordem com coeficientes constantes via método de Kansa.
Uma equação diferencial ordinária de segunda ordem com coeficientes
constantes é uma equação do tipo:
(4.36)
Onde é uma função de e são constantes reais.
Façamos para ilustrar . Além disso, tome
então temos, a seguinte equação diferencial ordinária:
(4.37)
Como condição de contorno, admitiremos que e Como
esta equação tem solução analítica e trata-se de uma EDO de 2ª ordem não
homogênea, podemos solucioná-la pelo método dos coeficientes a determinar e
compararemos sua solução analítica com a solução numérica, otimizando o
parâmetro de forma c. Lançando, para tanto, mão do método sem malhas
(MESHLESS). Então assim, primeiro passo é resolver a equação característica:
(4.38)
A função complementar é:
(4.39)
Uma solução particular para o problema é:
(4.40)
48
Então façamos:
(4.41)
Concluímos que
e , com isso a solução geral da EDO, é a soma da
função complementar com a solução particular:
(4.42)
(4.43)
Aplicando as condições de contorno e , segue-se que:
(4.44)
(4.45)
Logo a solução analítica é:
(4.46)
Modelaremos o problema da EDO de 2ª ordem em rotina do SCILAB. Para
tanto seguiremos alguns passos, apresentados no Fluxograma 4.2.1.
49
4.2.1 Fluxograma do processo de uso do método sem malhas
1. Identifica-se as constantes, coeficientes da equação
2. Delimita-se o intervalo relativo ao domínio da equação
3. Discretiza-se o domínio
4. Indica-se as condições de contorno do problema
5. Determina-se a matriz correspondente a , das distâncias
6. Determina-se a matriz FI
7. Obtém-se a matriz LFI
8. Monta-se a matriz de colocação A
9. O programa verifica se a matriz possui inversa
10. Determinam-se os valores de
11. Otimiza-se o parâmetro de forma
12. Indica-se o método: direto ou indireto
13. Plota-se os resultados em gráficos.
50
A rotina do problema em questão está apresentada no APÊNDICE B,
quando executada, calcula a solução aproximada do problema a partir da rotina de
cálculo desenvolvida no SCILAB. Após a execução da rotina em SCILAB solicita-se
ao usuário os seguintes dados por meio das caixas de diálogo a seguir:
Definição das condições de contorno e indicação do método: Número um,
direto (multiquadráticas) ou número dois método inverso (multiquadráticas inversas):
Figura 4.4 – Definição das condições de contorno da EDO.
Discretização do Domínio, temos que definir o intervalo considerado com o
limite inferior e o limite superior e numerarmos os pontos, para o refinamento dos
mesmos.
Figura 4.5 – Discretização do domínio da EDO.
51
Identificação das constantes ou coeficientes da equação diferencial, ou
seja, definição do tipo de EDO:
Figura 4.6 – Definição dos coeficientes da EDO.
A partir daí a rotina solicita a indicação do parâmetro de forma “c”, então
definimos um intervalo atribuindo um valor mínimo e máximo e um incremento e,
com isso, o programa busca um valor otimizado para o “c”, ou seja, um valor que
resultará na melhor aproximação possível entre a solução numérica e a solução
analítica do problema descrito.
Figura 4.7 – Dados de entrada para a otimização do parâmetro de forma c.
52
Figura 4.8 – Gráfico de resíduo do contorno e determinação de c = 68.21 como sendo o c - ótimo, via método multiquadrático direto.
Figura 4.9 – Gráfico de resíduo do domínio e determinação de c=15.81, como sendo o c - ótimo, via método multiquadrático direto.
53
Figura 4.10 – Comparação entre solução numérica e solução analítica, pontos do domínio, com c=15.81, otimizado pelo método multiquadrático direto.
Na indicação do número dois na caixa de diálogo abaixo e seguindo as
mesmas opções já descritas, encontraremos a aproximação pelo método inverso, da
solução numérica com a solução analítica:
Figura 4.11– Definição das condições de contorno da EDO, para uso do método inverso.
54
Figura 4.12 – Gráfico de resíduo e determinação de c=70.61, como sendo o c - ótimo do contorno, via método multiquadrático inverso.
Figura 4.13 – Gráfico de resíduo e determinação de c=16.81, como sendo o c - ótimo do domínio, via método multiquadrático inverso.
55
Figura 4.14 – Comparação entre solução numérica e analítica, pelo método inverso c=16.81.
4.3 Equação da difusão pura
A equação da difusão é descrita como,
(4.47)
Determinaremos a solução da equação da difusão, fazendo uso do método
de separação de variáveis e de série de Fourier, obedecendo às condições iniciais e
de contorno expressas, respectivamente, abaixo,
{
(4.48)
Sendo que D representa o coeficiente de dispersão. Façamos a separação de
variáveis, supondo:
(4.49)
56
Assim, deriva-se a função , em relação a ,
(4.50)
Determinamos a segunda derivada da função , em relação a ,
(4.51)
Separamos as variáveis, igualando posteriormente a uma constante – :
(4.52)
(4.53)
{
(4.54)
Resolvemos as equações diferenciais ordinárias:
{
(4.55)
{
(4.56)
Resolvendo a equação característica da primeira equação do sistema (4.56), tem-se
que:
(4.57)
57
√
√
(4.58)
Concluímos que:
(
√ ) (
√ )
(4.59)
Usando as condições de contorno,
(4.60)
(
√ )
(4.61)
(4.62)
(4.63)
∫
∫
(4.64)
(
)
(4.65)
Usa-se a série de Fourier,
∑
(
)
(4.66)
∫ (
)
(4.67)
58
Tomando, em particular,
(4.68)
(4.69)
(4.70)
∫
(4.71)
(4.72)
A solução analítica, abaixo, representa a concentração da
substância em função do espaço percorrido e do tempo.
∑
(4.73)
Usando a abordagem Kansa na equação da difusão, temos a seguinte
aproximação:
∑
(4.74)
Na forma compacta (matrizes), podemos escrever
(4.75)
Onde é definido por:
√
(4.76)
59
Nesse caso temos uma função unidimensional. Escolhendo o método direto
multiquadrático, , como função de base, temos:
√
√
(4.77)
Sendo que é um parâmetro de forma. De acordo com a equação da difusão será
necessário calcular as derivas parciais em relação a de primeira e segunda ordem.
Calculando as derivadas parciais temos:
•Derivada parcial de primeira ordem é igual a derivada ordinária de função
de uma variável,
√
(4.78)
• Consequentemente a derivada parcial de segunda ordem, também será
igual a derivada de segunda ordem ordinária de função de uma variável,
√( )
(4.79)
Sabendo que é função apenas de ou seja, e que é função
do tempo, ou seja, , substituindo na equação da difusão (equação 4.46)
ficamos com:
(4.80)
Usando um esquema explícito para o tempo, ficamos com:
(
)
(4.81)
Como o termo desconhecido é , isolando, ficamos com a equação:
60
(
)
(4.82)
Considera-se que:
(
)
(4.83)
Assim, temos a equação básica do método meshless para a equação da difusão
será:
(4.84)
Para o método sem malhas ficar completo, devemos introduzir as
condições iniciais e de contorno no nosso problema. Inserindo as matrizes para
, ficamos com:
(4.85)
A matriz é chamada de matriz de colocação. A matriz é a
concatenação da primeira linha da matriz e última linha da matriz com a matriz
, ou seja:
[
]
(4.86)
e será a concatenação das condições de contorno e condições inciais,
com isso, segue-se que:
[
] [
]
(4.87)
61
A partir dessa formulação podemos obter para e
consequentemente os outros e assim a solução numérica em qualquer intervalo de
tempo. A solução numérica em será dada por:
(4.88)
onde nos vimos como é feito o cálculo.
Os demais termos de serão obtidos a partir da seguinte equação:
(4.89)
Onde será dado por:
(4.90)
A estabilidade do sistema será mantida pela grandeza adimensional, conhecida
como número de difusão :
(4.91)
De acordo com Carnahan, (1969) para
assegura-se que a solução
não vai oscilar, e
tende a minimizar os erros de truncamento. Assim definindo
um , o passo de tempo será:
(4.92)
Vamos agora otimizar o valor do parâmetro de forma . A otimização de
foi feita com base nos resíduos, e foi considerado que não varia com o passar do
tempo. Para cada valor de , obtemos um vetor resíduo (para ). O vetor resíduo
é dado por:
(4.93)
62
Será mais conveniente trabalhar com o somatório do quadrado desse
resíduo. Logo para cada tem-se um somatório, aquele menor corresponde ao
melhor . Após ter encontrado o melhor , deve-se encontrar os outros , e
consequentemente a solução numérica em qualquer tempo. Conforme os gráficos,
abaixo:
Figura 4.15 – Solução analítica da Equação da Difusão Pura para os instantes t=0s; t=0.05s; t=0.10s com N=50 pontos, x (Km).
63
Figura 4.16 – Comparação solução analítica com solução numérica da Equação da Difusão Pura, pontos do domínio, para os instantes t=0s; t=0.05s; t=0.10s com N=5 pontos, x (Km).
Figura 4.17 – Comparação solução analítica com solução numérica da Equação da Difusão Pura, pontos do domínio, para os instantes t=0s; t=0.05s; t=0.10s com N=20 pontos, x (Km).
64
Figura 4.18 – Comparação solução analítica com solução numérica da Equação da Difusão Pura, pontos do domínio, c-ótimo 0.3, para os instantes t=0s; t=0.05s; t=0.10s com N=40 pontos.
Figura 4.19 – Gráfico comparativo de c e do somatório do quadrado do resíduo
65
4.4 Equação Diferencial Parcial da Difusão Advecção
A equação da difusão advecção é dada por:
(4.94)
Onde é a velocidade, representa a dispersão e é a concentração.
As condições iniciais são descritas pela função delta de Dirac e isto
significa que a concentração no instante inicial é zero em todos os pontos,
com exceção do segmento de reta , no qual foram despejados uma
quantidade de um determinado poluente com massa em gramas. A função delta
de Dirac rege o problema da massa por unidade de área :
(4.95)
As condições de contorno do problema podem ser pensadas de modo que
a partir de uma “longa distância” a concentração seja anulada, ou em termos mais
preciso, no limite temos:
(4.96)
Eis a solução analítica da equação de difusão advecção (FISHER, et al.1979, apud
ANDRADE, 2006):
√ (
)
(4.97)
66
4.4.1 Problema envolvendo a equação da difusão advecção
Um derramamento de uma substância conservativa, poluente, com massa
equivalente a ocorre em um rio hipotético a uma distância a jusante de
uma cidade. A profundidade do rio é de 1m, a largura é de , a velocidade
é de e a dispersão equivale a . Representar graficamente,
a distribuição da concentração da substância conservativa ao longo do suposto rio,
em um tempo de fazendo a comparação entre a solução numérica, via
método sem malhas e a solução analítica apresentada. Em seguida, compare em
um mesmo gráfico a solução numérica e a solução analítica da distribuição da
concentração da substância, nos seguintes instantes:
a) 12 minutos;
b) 24 minutos;
c) 36 minutos.
Onde é massa total das partículas normais a área da seção transversal, ou seja,
e representa a distância entre a cidade e o ponto de lançamento.
(4.98)
(4.99)
Em cada caso é dado por:
(4.100)
67
4.4.2 Abordagem Kansa na equação de difusão advecção
Usando a abordagem Kansa na equação de Difusão Advecção, temos a
seguinte aproximação:
∑
(4.101)
Onde é definido por:
√
(4.102)
Nesse caso, temos uma função unidimensional. Escolhendo como função de
base, temos:
√
√
(4.103)
De acordo com a equação da difusão será necessário calcular as derivas parciais
em relação a de primeira e segunda ordem. Calculando as derivadas temos:
• Derivada primeira
√
(4.104)
• Derivada segunda
√( )
(4.105)
68
Para uma variação temporal, um esquema implícito pode ser utilizado.
Aplicando a forma implícita na equação da difusão teremos:
(4.106)
onde indica o intervalo de tempo e o subescrito é o valor desconhecido (ou
novo), a ser resolvido. A solução aproximada pode ser expressa por:
∑
(4.107)
Os outros termos da equação (4.106) são:
∑
(4.108)
∑
(4.109)
∑
(4.110)
Substituindo esses termos na equação (4.106) obtêm-se:
∑
∑
∑
∑
(4.111)
Isolando o termo desconhecido , a equação fica:
∑
∑
∑
∑
(4.112)
69
Multiplicando ambos os membros por e pondo em evidencia no membro
esquerdo da equação, segue-se que:
∑
*
+
∑
(4.113)
Para , que produz um sistema linear de equações para o valor
desconhecido .
As matrizes que representam a equação (4.110) são:
(4.114)
[ (
) (
)
(
) (
)
(
) (
)
]
(4.115)
[
]
[
]
[
]
Dessa forma podemos reescrever a equação (4.110) em sua forma matricial, que é
dada por:
(4.116)
Onde , , e são matrizes. Como queremos encontrar a matriz ,
devemos isolar esse termo, para isso, multiplicamos pela inversa da matriz , ou
seja:
70
(4.117)
Da Álgebra de matrizes, sabemos que uma matriz multiplicada pela sua inversa é
igual à matriz identidade.
(4.118)
Sendo que é a matriz identidade N-dimensional definida por:
[
]
(4.119)
Fazendo algumas operações algébricas, obtemos:
(4.120)
4.4.3 Matrizes da equação de difusão advecção
Estando a forma geral da equação pronta, o próximo passo é montar as
matrizes que caracteriza o método sem malhas e, em seguida, implementaremos na
plataforma SCILAB a rotina para que a mesma seja executada . As matrizes são:
matriz , matriz , matriz
e matriz
. De acordo com a equação (4.106) a matriz
é uma composição destas matrizes. As matrizes são:
• Matriz
[√
] [| |] (4.121)
71
[
]
[
]
(4.122)
• Matriz
[√
] [√ ]
(4.123)
[ √
√ √
√ √
√
√ √
√
]
(4.124)
• Matriz
[
√
]
*
+
(4.125)
[
]
[
]
(4.126)
72
• Matriz
[
(√ )
] [
] (4.127)
[
]
[
]
(4.128)
Com as matrizes principais montadas, podemos agora determinar a matriz
de colocação . Cada termo dessa matriz será dado por:
(4.129)
Na forma matricial temos:
[ (
) (
)
(
) (
)
(
) (
)
]
(4.130)
73
Substituindo os valores das matrizes principais teremos: (4.131)
[ (
) (
) (
)
(
) (
) (
)
(
) (
) (
)
]
Onde e são constantes e é o intervalo de tempo.
Logo abaixo, figura 4.20, está exposto o fluxograma, simplificado, do
processo de comparação entre solução numérica e solução analítica com otimização
do “c” no scilab:
Figura 4.20 – Fluxograma simplificado da otimização do c.
74
4.4.4 Estabilidade do sistema Número de Courant
Para Garantir a estabilidade do sistema foi adotado o número
adimensional Número de Courant. Essa grandeza adimensional representa, então,
por meio de uma razão de distâncias, o fluxo advectivo que atravessa os volumes da
malha em um dado intervalo de tempo. Se supormos uma malha e uma velocidade
uniformes numa região do domínio, o Número de Courant local torna-se o número
de volumes da malha atravessados por uma perturbação naquele passo de tempo é
dado por Carnahan (1996):
(4.132)
Onde é a velocidade é o passo de tempo e é o passo no espaço. Para um
sistema estável esse número deve ser tal que:
(4.133)
para a equação da difusão foi adotado que , dessa maneira a partir de um
passo podemos encontrar o passo de tempo. Isolando e usando o
passo de tempo utilizado foi:
(4.134)
75
Figura 4.21 – Gráfico comparativo da solução numérica e da solução analítica da Equação de Difusão-Advecção com t=15 min e N=100 pontos.
Figura 4.22 – Gráfico comparativo da solução numérica e da solução analítica da Equação de Difusão-Advecção com t=15 min e N=200 pontos.
76
Figura 4.23 – Gráfico comparativo entre a solução numérica e a solução analítica da Equação de Difusão-Advecção, pontos do domínio, nos instantes t=12min, t=24min, t=36min. O parâmetro de forma, c = 0.149254, foi otimizado pelo SCILAB 5.4.1.
Figura 4.24 – Gráfico comparativo do parâmetro de forma e do resíduo da equação de Difusão-Advecção. Para cada c, tem-se um resíduo e o c-ótimo implica em um menor resíduo.
77
5. CONCLUSÕES:
Nesse trabalho foi dada uma atenção às equações diferenciais, ordinárias
e parciais, haja vista uma vasta gama de problemas físicos e, por conseguinte de
engenharia ser regidos por estas equações. Assim foram aplicados métodos sem
malha, funções de base radial, multiquadráticas e multiquadráticas inversas com
aproximação de Kansa e, com isso, foram feitas simulações, utilizando a plataforma
SCILAB e comparando, graficamente, a aproximação entre a solução numérica e a
solução analítica das supracitadas equações.
Nas simulações utilizamos, a priori, uma equação diferencial linear de
primeira ordem, com condição inicial, para ilustrar o método sem malhas e comparar
a solução numérica com a solução analítica. Nesse primeiro caso, nós usamos um
parâmetro de forma , suposto conforme literatura (FERREIRA, p.87) como a razão
de dois pela raiz quadrada de , onde é o número de pontos.
Em seguida, apresentamos um problema de valor de contorno em uma
equação diferencial linear de segunda ordem, não homogênea, onde o parâmetro de
forma , foi otimizado via função de base radial multiquadráticas e multiquadráticas
inversa, por meio da plataforma SCILAB.
Comparamos a solução numérica com a solução analítica da equação
diferencial parcial da difusão pura e, com isso, encontrando um parâmetro de forma
c ótimo, minimizando o erro e, consequentemente o resíduo da aproximação.
Ao final, simulamos um problema de contaminação de um determinado
aquífero por meio de uma substância conservativa. A modelagem matemática do
problema foi descrita por meio de uma equação “1D transiente,” denominada de
equação de difusão-advecção, onde comparamos a solução numérica com a
solução analítica e plotamos os resultados graficamente para analisarmos a
concentração do poluente em função do tempo e da distância percorridos.
A modelagem matemática, em conjunto com a simulação computacional e
o advento dos métodos numéricos, principalmente o método numérico sem malhas,
têm demonstrado ser eficientes e eficazes, nos estudos de recursos hídricos. Dessa
78
forma, são utilizados para prever situações nas quais estudos empíricos, sejam
dispendiosos e, por vezes também onerosos.
A vantagem do método sem malhas, frente aos demais métodos de
aproximação de uma equação diferencial está basicamente na relativa facilidade
com a qual podemos montar um algoritmo, pois não há necessidade de construção
de uma malha regular, uma vez que somente nos interessa a distância entre os
pontos discretizados e não a sua posição geométrica. Ademais, com o aumento da
dimensão, não há um incremento considerável, como ocorre no método dos
elementos finitos, na dificuldade de modelarmos problemas, usando como
ferramenta o SCILAB. Contudo, por se tratar de um ramo emergente vale pesquisar
maneiras de diminuir o esforço computacional, o qual ainda é grande, por conta do
número de operações que a Álgebra de matrizes e sistemas lineares exige.
Por fim, para estudos futuros podemos modelar diversos problemas de
Recursos Hídricos, tais contaminação de aquíferos, recarga hidráulica e golpe de
aríete, fazendo uso dos métodos sem malha, das funções de base radial, mas não
se limitando as funções multiquadráticas diretas e inversas, como fizemos nesse
trabalho e, sim fazendo uso da Gaussiana.
79
REFERÊNCIAS BIBLIOGRÁFICAS
AZEVEDO NETTO, J. M.; FERNANDEZ, M. F.; ARAÚJO, R.; ITO, A E. Manual de Hidráulica 8a ed. – São Paulo – 1998. 670 p.
BAPTISTA, Márcio Benedito; COELHO, Márcia Maria Lara Pinto. Fundamentos de Engenharia Hidráulica 2a ed. rev. – Belo Horizonte – Editora: UFMG, 2003. 440 p.
Barret R., Berry M., et. al. Templares for the solution of Linear Systems: Building Blocks for Iterative Methods. SIAM, 1994.
BENNET, C. O. & MYERS, T. E. Fenômenos de Transporte. São Paulo, McGraw-Hill do Brasil.1980.
BOYCE, William Edward. Equações diferenciais elementares e problemas de valores de contorno.9.ed. São Paulo: LTC, 2010.
BRAGA, Benedito, et.al. Introdução à Engenharia Ambiental. São Paulo, SP: Pearson Prentice Hall, 2005.
BRANCO, Samuel Murgel, Água: Origem, uso e preservação.2.ed. São Paulo:
Moderna, 2003.
BRANNAN, J. R.; BOYCE, W. E. Equações Diferenciais: uma introdução aos métodos modernos e suas aplicações. 1ª ed. Rio de Janeiro: LTC, 2009.
CARNAHAN, B., LUTHER, H. A. e WILKES, J. O. (1969) Applied Numerical Methods. John Wiley & Sons, Inc.
CHOW, Ven Te; MAIDMENT, R. David; MAYS, Larry W.; Hidrologia Aplicada. Tradução para o espanhol de SALDARRIAGA, Juan G.; Santa Fé de Bogotá- Colômbia, 1994. Editora: Mc GRAW-HILL.
COURANT, R. Variational methods for the solution of problems of equilibrium and vibrations, Bulletin of American Mathematical Society, 1943.
DAKER, Alberto. Captação, elevação e melhoramento da água; 7a ed. Volume II – Rio de Janeiro - 1987. 408 p.
DIAS, Nelson L; Obtenção de uma solução analítica da equação de difusão-
advecção com decaimento de 1ª ordem pelo método da transformação de similaridade generalizada. RBRH, 2002, vol. 8 n.1 p.181-188
DUARTE, C. Armando, A Review of Some Meshless Methods to Solve Partial Differential Equations, Tech Report, TICAM – Texas Institute for Computational and Applied Mathematics.1995.
FERREIRA, A. J. M.; KANSA, E.J et. al. Progress on Meshless Methods. Porto. Portugal: Springer 2009.
FISCHER, H. B., et. al.: Mixing in Inland and Coastal Waters. Academic Press, New York 1979 apud ANDRADE, Rodrigo Campos: Uma nova abordagem para a
solução numérica de problemas de advecção e difusão multidimensional em corpos de água naturais. Tese de doutorado. UFRJ/ COPPE, 2006
80
FOSTER, S. S. D. ; HIRATA, R. C.; ROCHA, G.A. 1998. Riscos de poluição de
águas subterrâneas: uma proposta metodológica de avaliação regional. In: CONGRESSO BRASILEIRO DE ÁGUAS SUBTERRÂNEAS., 5. São Paulo. Anais. São Paulo, ABAS.
GADELHA FILHO, D. Apostila SCILAB 5.1.1., Universidade Federal do Ceará, 2004
GARCEZ, Lucas Nogueira. Elementos de Mecânica dos Fluidos: Hidráulica Geral. Ed. Edgard Blucher, 1963.
GILES, Ronald V. Mecânica dos Fluídos e Hidráulica. São Paulo, McGraw-Hill do Brasil Ltda.1982.
GILES, Ronald V. Mecânica dos Fluídos e Hidráulica. São Paulo, McGraw-Hill do Brasil Ltda, 1996.
IÓRIO, Valéria. EDP, um curso de graduação. 3 ed. Rio de Janeiro: IMPA, 2010.
LIMA, Elon Lages. Álgebra Linear. 8.ed.Rio de Janeiro: Impa, 2011. (Coleção matemática universitária).
LIU, Gui-Rong. Mesh free methods: moving beyond the finite elemento method. USA, New York: CRC Press LLC, 2003.
MACHADO, Kleber D. Equações Diferenciais Aplicadas à Física. 2ª ed. Ponta Grossa, EDUEPG, 2003.
MENESCAL, Germana Cavalcante. Modelagem, Numérico – Analítica da contaminação de aquíferos utilizando o método de colocação RBF livre de malha. 128f. Tese (Doutorado em Engenharia Civil – Recursos Hídricos) –
Departamento de Engenharia Hidráulica e Engenharia Ambiental, Universidade Federal do Ceará – UFC, Fortaleza, 2008.
MEYER, J.F.C.A.; SALVATIERRA, M.M; Modelagem matemática e simulação computacional da presença de materiais impactantes tóxicos em casos de
dinâmica populacional com competição inter e intra-específica. Revista Biomatemática, 2006, vol.16, p. 1-22
MIJARES, Francisco Javier Aparicio; Fundamentos de hidrologia de superfície. México: Limusa, 1992.
MUNSON, Bruce R.; DONALD F. Young; THEODORE, H. Okiishi. Fundamentos de
mecânica dos Fluidos. São Paulo: Edgar Blücher, 2004.
NETTO, J. M. de A. Manual de Hidráulica. São Paulo, McGraw-Hill do Brasil Ltda, 1990.
NEVES, Eurico Trindade. Curso de Hidráulica, 9a ed.- São Paulo – 1989. 577 p.
PHILIPPI, Arlindo Jr; ROMÉRO, Marcelo de Andrade; BRUNA, Gilda Collet. Curso
de Gestão Ambiental. Barueri, SP: Manole, 2004.
PIRES, Paulo Sérgio da Motta, Apostila de Introdução ao Scilab - Versão 3.0. Departamento de Engenharia de Computação e Automação, Universidade Federal
81
do Rio Grande do Norte-UFRN. Natal-RN, Julho de 2004.
http://www.dca.ufrn.br/~pmotta. PORTO, Rodrigo de Melo, Hidráulica básica.4.ed. São Carlos: EESC-USP,2006
SHAMES, Irving H. Mecânica dos Fluídos, vols. 1 e 2, São Paulo, Edgard Blücher.1980.
STEINBRUCH, Alfredo; WINTERLE, Paulo. Álgebra linear: teoria e problemas. 2.ed. São Paulo: Makron, 1987.
STREETER, Victor L. Mecânica dos Fluídos. São Paulo, McGraw-Hill do Brasil Ltda.1980.
VIANNA, Marcos Rocha; Mecânica dos Fluidos para Engenheiros - 4a ed. Belo Horizonte – 2001. 582 p.
ZILL, Dennis G. and CULLEN, Michael R., Differential Equations with Boundary-
Value Problems; Seventh Edition, Pearson, New York, 2007;
ZILL, Dennis G. e Cullen, Michael R. Equações Diferenciais, 3ª Edição, Volume 1, Pearson S.A, 2001.
ZILL, Dennis G. Equações diferenciais com aplicações em modelagem. São
Paulo: Pearson, 2011.
82
APÊNDICE A – NOCÕES DE DEPENDÊNCIA LINEAR E INDEPENDÊNCIA
LINEAR E FUNÇÃO DE IMPULSO DE DIRAC.
Uma propriedade importante de conjuntos e que influencia à formação de
bases de espaços vetoriais é a identificação da dependência e da independência
lineares. Sabe-se através da literatura (LIMA, 2011) que um conjunto de vetores,
forma uma base de um Espaço Vetorial, quando os mesmos são linearmente
independentes e geram o espaço. Nesse sentindo, existe uma ampla gama de
espaços vetoriais, como o , como o espaço das matrizes quadradas de ordem n,
como o espaço de polinômios de grau menor ou igual a n, espaço dos operadores
lineares de dimensão finita, entre outros.
A.1 – Conjunto linearmente dependente (LD)
Seja V um espaço vetorial sobre o corpo , e sejam vetores de
V. o conjunto , é dito linearmente dependente sobre se existem
escalares reais ; não todos nulos, tais que:
(A.1.1)
Para todo .
Os vetores e são LD, desde que .
Digamos , com então temos, conforme figura A.1.1, uma expansão do
vetor .
Figura A.1 – Expansão do vetor u
83
Se , temos uma mudança de sentindo, digamos , então
representamos, geometricamente, através da figura A.2:
Figura A. 2 – Mudança de sentido do vetor u
Caso o | | , haverá uma contração que se dará no sentido do vetor
se , e no sentindo inverso se Em qualquer caso, quando um vetor é
múltiplo do outro, dizemos que os mesmos são LD.
A.2 – Conjunto linearmente independente (LI)
Seja V um espaço vetorial sobre o corpo , e sejam vetores de
V. o conjunto , é dito linearmente independente sobre se existem
escalares reais ; tais que:
(A.2.1)
Se, e somente se, =0, para todo .
Como ilustração podemos citar a base canônica do
(A.2.2)
84
A combinação linear, abaixo resulta que =0
para todo :
(A.2.3)
A.3 – Aplicação Linear (Transformação Linear)
Sejam dois espaços vetoriais e , sobre o corpo dos reais, . Uma
aplicação linear é uma função , que obedece às seguintes propriedades:
i) ( ( ; (A.3.1)
ii) ( ( , . (A.3.2)
A.4 – Aplicação linear do vetor nulo entre espaços vetoriais.
Sendo A uma aplicação linear de dois espaços vetoriais e , sobre o
corpo dos reais, . A aplicação Linear:
( ) (A.4.1)
De fato, note que como é espaço vetorial, temos que , então
( ) . Por conseguinte, ( ) ( ), mas como A é uma aplicação linear,
segue-se que: ( ) ( ) ( ) ( ) Conclui-se, assim, a
demonstração.
85
A.5 – Funçõesde Impulso
Alguns fenômenos tem natureza impulsiva como, por exemplo, os
problemas de voltagem, ou forças muito grandes em intervalos de tempo muito
curto. Esses problemas geralmente são modelados através de equações diferenciais
da forma:
(A.5.1)
Onde é grande em um intervalo pequeno e é zero nos outros
pontos. A integral, definida abaixo por:
∫
(A.5.2)
Ou, como fora do intervalo de tempo ,
∫
(A.5.3)
Funciona como uma medida de força do termo não homogêneo, ou seja
em sistemas mecânicos, onde a função é uma força, é o impulso total da
força g(t) sobre o intervalo no qual a força está definida . Em
condições similares, no caso de um circuito elétrico e é a derivada em relação
ao tempo da voltagem, então a integral representa a voltagem total no circuito
no intervalo aberto considerado .
A função impulso unitário denominada de função delta ( ) de Dirac tem as seguintes
propriedades:
i. sempre que ;
ii. ∫
Vale ressaltar que a função delta de Dirac, na verdade não é uma função, no
sentindo clássico ou formal da palavra, mas é conhecida como à classe de
funções generalizadas.
86
Suponha, em particular que ,
Sendo que é uma constante positiva pequena.
Figura A.3 – Gráfico da função delta de Dirac de .
Fonte: Boyce (2010 p. 184)
Figura A.4 – Gráfico da função delta de Dirac de , quando .
Fonte: Boyce (2010 p. 184)
87
APÊNDICE B – Breve exposição sobre o SCILAB
Figura B. 1– Sobre o SCILAB 5.4.1
B.1 – INTRODUÇÃO
O SCILAB é um freeware matemático, ou software livre, que faz um
grande número de operações matemáticas e é adequado para a resolução de
problemas numéricos de matemática e, por conseguinte de Física e Engenharia.
Com isso, proporciona um ambiente de computação muito frutífero para aplicações
em uma ampla gama de problemas de Engenharia Civil, em particular na
Engenharia de Recursos-Hídricos.
O SCILAB é uma versão livre com mesma plataforma do MATLAB,
software muito utilizado por cientistas, matemáticos e engenheiros. Este programa
foi lançado (PIRES, 2004) em código aberto, sendo criado e mantido por
pesquisadores do INRIA (Institut de Recherche em Informatique et em Automatique),
através do Projeto MÉTALAU (Méthods, algorithmes et logiciels pour l’automatique)
e também com o respaldo da ENPC ( École Nationale des Ponts et Chaussées).
O SCILAB pode ser instalado nos seguintes sistemas operacionais
GNU/Linux, Mac OS X e Windows XP/Vista/7/8, reconhece vários tipos de
linguagens e códigos como o “c”, ou “c++” e o “FORTRAN”, entre outros. Além disso,
gera gráficos 2D e 3D e executa simulações em sistemas mecânicos, hidráulicos e
88
elétricos onde estão localizados os arquivos executáveis. Nas distribuições
Unix/Linux o script SCILAB aciona o executável SCILEX. Neste diretório estão,
também, localizados programas para manipulação de arquivos POSTSCRIPT e
LATEX gerados pelo SCILAB.
Entre as várias funcionalidades desse software, estão as simulações
matemáticas, visualização de gráficos, optimização, estatística, álgebra linear,
cálculo diferencial e integral, EDO e EDP. A versão utilizada nesse trabalho foi a
scilab-5.4.1, para Windows 8, PC 64 bits.
Ao abrirmos o SCILAB nos deparamos com o console, ambiente virtual
onde são inseridos os comandos, preferencialmente curtos para que sejam logo
executados.
Figura B.2 – Ambiente, plataforma SCILAB 5.4.1
Figura B.3 – Acessando o aplicativo Scinotes do SCILAB 5.4.1
89
Programações mais extensas, como as desenvolvidas nesse trabalho,
devem ser feitas no ambiente virtual Scinotes, para a posteriori serem executadas.
Figura B.4 – Aplicativo SciNotes do SCILAB 5.4.1
B 1.1 – Criação de Matrizes
Uma matriz é uma espécie de tabela que permite guardar números de uma
forma ordenada e indexável, a mesma também pode ser vista como vetor, haja vista
o conjunto de matrizes m por n, satisfazer os oito axiomas de espaço vetorial. No
SCILAB para criarmos uma matriz 3x3, por exemplo, basta digitá-la, separada por
ponto e vírgula e colocada entre colchetes, conforme feito na FIGURA B.5:
Figura B.5 – Gerando uma matriz 3x3 no SCILAB 5.4.1
90
B 1.2 – Dimensão de Matrizes
A dimensão da matriz A é determinada no SCILAB, pelo comando size(A),
à qual nesse caso, em particular, conforme figura abaixo, é três por três.
Figura B.6 – Determinando a dimensão de uma matriz A, size(A), no SCILAB 5.4.1
91
B 1.3 – Transposição de Matrizes
Transpor uma matriz é transformar linhas em colunas e,
consequentemente, colunas em linhas. Então assim, a transposta da matriz A,
apresentada na FIGURA B.7 é dada por C=A’.
Figura B. 7 – Cálculo da transposta de uma matriz A, size(A), no SCILAB 5.4.1
92
B 1.4 Inversa de Matrizes
O cálculo da inversa de uma matriz invertível, ou seja, cujo determinante é
diferente de zero, utilizando o SCILAB é muito simplificado, conforme demonstramos
abaixo, bastando apenas, usarmos a representação inv(A) para representarmos a
matriz inversa de A.
Figura B.8 – Cálculo da inversa de uma matriz A, no SCILAB 5.4.1
93
B 1.5 – Comandos para Iterações (for e while)
Existem dois comandos que permitem a realização de iterações, loops,
implementados através do comando for, ou fazendo uso do comando while. A
instrução, usando o comando for está representada na FIGURA B.9, abaixo:
Figura B. 9 – Esquema de uso do comando “for” no SCILAB 5.4.1
94
Figura B.10 – Aplicação do comando “for” no SCILAB 5.4.1
O comando loop baseado no while realiza uma sequência de instruções
enquanto uma determinada condição tiver que ser satisfeita, geralmente está incluso
à comparação de objetos, tais como os operadores descritos na tabela abaixo:
Fonte: (PIRES, 2004, p. 59)
Tabela B.1 - Operadores básicos do SCILAB
Operadores Significado
igual a
menor do que
maior do que
menor ou igual
maior ou igual
diferente
95
Figura B.11 – Aplicação do comando “while” no SCILAB 5.4.1
96
B 1.6 – Comandos e funções pré-definidos no SCILAB
O SCILAB é um programa versátil, haja vista a imensa quantidade de
problemas, os quais podem ser trabalhados e, portanto, modelados com essa
ferramenta. Ademais, o programa oferece um recurso de ajuda, no qual, pode-se ter
acesso a funções pré-definidas em diversos ramos da matemática como: Álgebra
Linear, Equações diferenciais, polinômios, cálculo diferencial, estatística, entre
outros.
Figura B.12 – Navegador de ajuda no SCILAB 5.4.1, possui comandos específicos para os mais variados ramos da Matemática.
97
APÊNDICE C - ROTINA DO PROBLEMA 4.1 DA EDO DE 1ª ORDEM.
//Dados de entrada--------------------------------------------------------------
N=evstr(x_dialog('Insira o número de pontos (deve ser um inteiro): ','7')) X1=0; XN=2;
Y1=4 C=2/(N)^0.5; //------------------------------------------------------------------------------
//Criação de matrizes e vetores nulos-------------------------------------------
X=[];
XI=[]; XJ=[]; R=[];
FI=[]; LFI=[]; A=[];
LAMBDA=[]; FF=[]; solA=[];
solN=[]; //------------------------------------------------------------------------------
//Intervalo---------------------------------------------------------------------
intervalo=(XN-X1)/(N-1) for (i=1:1:N)
X(i)=intervalo*(i-1)+X1 end //------------------------------------------------------------------------------
//XI(Linha) --------------------------------------------------------------------
for (linha=1:1:N)
XI(linha,1)=X(linha) end //------------------------------------------------------------------------------
//XJ(Coluna)--------------------------------------------------------------------
for (coluna=1:1:N)
XJ(1,coluna)=X(coluna) end //------------------------------------------------------------------------------
//Matriz R---------------------------------------------------------------------- for (i=1:1:N) //i para linha
for (j=1:1:N) //j para coluna
R(i,j)=(XI(i,1)-XJ(1,j)) end
end
98
//------------------------------------------------------------------------------
//Matriz FI---------------------------------------------------------------------
for (i=1:1:N)
for (j=1:1:N) FI(i,j)=(R(i,j)^2+C^2)^0.5 end
end //------------------------------------------------------------------------------
//Matriz LFI--------------------------------------------------------------------
for (i=1:1:N) for (j=1:1:N)
LFI(i,j)=(R(i,j)^2+C^2)^0.5+(R(i,j)/(R(i,j)^2+C^2)^0.5) end end //------------------------------------------------------------------------------
//Matriz A----------------------------------------------------------------------
for (i=1:1:N) for (j=1:1:N) if(i==1) then
A(i,j)=FI(i,j) else A(i,j)=LFI(i,j)
end end end //------------------------------------------------------------------------------
//FF----------------------------------------------------------------------------
FF(1,1)=Y1
for(i=2:1:N) FF(i,1)=X(i) end //------------------------------------------------------------------------------
//Inversão da matriz A e vetor LAMBDA------------------------------------------
invA=inv(A) LAMBDA=invA*FF //------------------------------------------------------------------------------
//Solução analítica e numérica--------------------------------------------------
solN=FI*LAMBDA
for (i=1:1:N) solA(i,1)=X(i)-1+5*exp(-X(i)) end //------------------------------------------------------------------------------
//Exibe os dados----------------------------------------------------------------
99
caminho=get_absolute_file_path('Meshless - Exemplo 1.sce')
files=caminho+'exibeDados.sci' exec(files, -1) //executa o programa no mesmo diretório exibe=%F //Exibe=falso
exibe=exibeDados(N,X,R,FI,LFI,A,solA,solN,invA,LAMBDA,FF); //se exibir os dados retorna verdadeiro (exibe=verdadeiro) //------------------------------------------------------------------------------
//Plotar gráfico formatado------------------------------------------------------
caminho=get_absolute_file_path('Meshless - Exemplo 1.sce')
files=caminho+'exibeGrafico.sci' exec(files, -1) //executa o programa no mesmo diretório exibeG=%F //Exibe=falso
estilo=[2,4] exibeG=exibeGrafico(estilo,X,solA,solN) //se exibir o gráfico retorna verdadeiro (exibe=verdadeiro) //------------------------------------------------------------------------------
100
APÊNDICE D - ROTINA DO PROBLEMA DA EDO DE 2ª ORDEM NÃO
HOMOGÊNEA.
//ROTINA PARA SOLUÇÃO NUMÉRICA DE EDO's UTILIZANDOO
//MÉTODO SEM MALHAS(MESHLESS) //EQ.GERAL:
//A*(d^2/du(x))+B*(d/du(x))+C*u(x)=D*sin(E*x)+F*cos(G*x)+H*x^3+I*x^2+J*x+L
//------------OBTENDO O PARÂMETRO DE FORMA "C" E DADOS PARA DISCRETIZAR O
DOMÍNIO
clear; //Apaga todas as variáveis da memória //++++++++++++++++++++++++++
//----------------------------------------
txt0=['Digite 1 para u(x1) ou 2 para Dx(x1)'; 'Condição de contorno 1 (cc1):';
'Digite 1 para u(x2) ou 2 para Dx(x2):';
'Condição Contorno 2 (cc2):'; 'Método direto 1, inverso 2'];
valor0=x_mdialog('Definindo Condições de Contorno',txt0,['1';
'0.0'; '1';
'0.0';
'1']); //Sugere valores
Escolhacc1=evstr(valor0(1));
cc1=evstr(valor0(2));
Escolhacc2=evstr(valor0(3));
cc2=evstr(valor0(4)); Metodo=evstr(valor0(5));
//------------------------------------------ //------------------------------------------
txt=['Número de pontos do domínio:';
'Limite inferior do domínio(x1):';
'Limite superior do domínio(x2):'; 'Parâmetro de forma C:']; //Indica nome das variáveis
valor=x_mdialog('Discretização do Domínio ',txt,['6.0';
'0.0'; '1.0';
'0.0']); //Sugere valores
np=evstr(valor(1)); xmenor=evstr(valor(2));
xmaior=evstr(valor(3));
PFindicado=evstr(valor(4));
101
if np<=1 then np=3//trata o caso do usuário colocar um único ponto
end
//------------DISCRETIZAÇÃO DO DOMÍNIO
deltaP=abs(xmaior-xmenor)/(np-1);//Calcula o incremento para discretização
Pd=xmenor:deltaP:xmaior; //Vetor linha com os pontos do domínio My=Pd';Mx=My; //Transforma vetor linha em etor coluna, armazena em Mx e My
//-----------OBTENDO VALORES DAS CONSTANTES DA EQUAÇÃO DIFERENCIAL
txt1=['A=';'B=';'C=';'D=';'E=';'F=';'G=';'H=';'I=';'J=';'L='];//Constantes EDO valor=x_mdialog('EDO: A.Dxx(u) + B.Dx(u) + C.u = D.sen(E.x) + F.cos(G.x) + H.x^3 + I.x^2 +
J.x + L',txt1,['1';
'0';'-9';'0';'0';'0';'0';'0';'0';'3';'0']); //Sugere valores
A=evstr(valor(1)); B=evstr(valor(2));
C=evstr(valor(3));
D=evstr(valor(4)); E=evstr(valor(5));
F=evstr(valor(6));
G=evstr(valor(7));
H=evstr(valor(8)); I=evstr(valor(9));
J=evstr(valor(10));
L=evstr(valor(11)); //(44-54)Armazena os valores das variáveis
//----------VERIFICA SE A EQUAÇÃO É DIFERENCIAL PARA INICIAR O MÉTODO
MESHLESS
if (A==0 & B==0) then //Plota um gráfico da função não diferencial plot(Mx,(D*sin(E*Mx)+F*cos(G*Mx)+H*Mx^3+I*Mx^2+J*Mx+L)*(1/C));
xlabel(["Pontos do Domínio - x"]);
ylabel("Valor da Função - y"); title("Gráfico da Função Não Diferencial");
legend(["Solução Exata"]);
disp('A equação não é diferencial.');//Avisa que a equação não é diferencial
abort; end
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++ //------BUSCA O VALOR ÓTIMO DE C PARA EDO NÃO-HOMOGÊNIA----COMEÇA AQUI
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++
//---------------OBTEM INFORMAÇÕES PARA OTIMIZAÇÃO DO PARÂMETRO FORMA if PFindicado==0.0 then
txt3=['C Inicial:';
'C Final:'; 'Incremento de C:']; //Indica nome das variáveis
valor=x_mdialog('Dados para otimização do C:',txt3,['0.01';
'100';
102
'0.20']); //Sugere valores
PFinicial=evstr(valor(1)); PFfinal=evstr(valor(2));
DeltaC=evstr(valor(3));
PF=PFinicial //---------INICIALIZA VARIÁVEIS
MPF=[] //Cria uma matriz vazia para armazenar os parâmetro C testado
MSR=[] //Cria uma matriz vazia para armazenar o resíduo testado
MRC= [] // Cria uma matriz vazia para armazenar os residuos das condições de contorno
k=0 //inicializa o contador
while(PF<PFfinal)
k=k+1 PF=PF+DeltaC
//------------DETERMINA A MATRIZ R [m n]=size(Mx) //obtem a quantidade de linha (m) e colunas (n)
for i=1:m //Laço para linhas
for j=1:m //Laço para as colunas
R(i,j)=([Mx(i,1)-My(j,1)]^2)^0.5; //Determina a matriz R end
end
//------------DETERMINA A MATRIZ FI [m n]=size(R) //obtem a quantidade de linha (ri) e colunas (rj)
for i=1:m //Laço para linhas
for j=1:n //Laço para as colunas
if Metodo == 2 then FI(i,j)=(R(i,j)^2+PF^2)^(-0.5);
else
FI(i,j)=(R(i,j)^2+PF^2)^0.5; //Determina a matriz FI end//end if
end //end for j
end// end for i
//------------CALCULA A MATRIZ L(FI)
[m n]=size(FI);
//Determina L(FI) para 2ª derivada for i=1:m //Laço para linhas
for j=1:n //Laço para as colunas
if Metodo == 2 then
RFI2(i,j)=(R(i,j)^2+PF^2)^(-2.5); LFI2(i,j)=RFI2(i,j)*((-(R(i,j)^2+PF^2))+(3*R(i,j)^2));//Matriz L(FI)
else
RFI2(i,j)=(R(i,j)^2+PF^2)^1.5; LFI2(i,j)=PF^2/RFI2(i,j); //Matriz L(FI)
end //end if
end//end for j
end//end for i
103
//Determina L(FI) para 1ª derivada
for i=1:m //Laço para linhas for j=1:n //Laço para as colunas
if Metodo == 2 then
LFI1(i,j)=-R(i,j)/(R(i,j)^2+PF^2)^(1.5);
else LFI1(i,j)=R(i,j)/(R(i,j)^2+PF^2)^0.5;
end//end if
end//end for j
end//end for i
//-------------MONTA A MATRIZ L(FI)
LFI=(A*LFI2)+(B*LFI1)+(C*FI); //Determina a matriz LFI
//-------------DETERMINA A MATRIZ 'A' E SUA INVERSA
if Escolhacc1 == 2 then
A1=LFI1(1,1:n); else
A1=FI(1,1:n);//recorta a 1ª linha da matriz FI (cc1)
end
A2=LFI(2:m-1,1:n); //recorta na matriz LFI (pontos interior do domínio)
if Escolhacc2 == 2 then//avalia se deve usar u ou du para cc2
A3=LFI1(m,1:n); //recorta a LFI1 (cc 2 "A derivada de FI")
else A3=FI(m,1:n);//recorta a última linha da FI (cc 2 u(x2)=u2.)
end
MA=[A1; A2; A3]; //Concatena linhas para formar a matriz A if det(MA)<>0 then
MAinv=inv(MA);
else disp('A matriz não possui inversa') //Avisa se a matriz não for invertível
abort
end
//-----------DETERMINA O VETOR Y
[m n]=size(Mx)
for i=1:m //Cria um laço descosiderando as condições de contorno da fronteira - MODIFIQUEI
AQUI y(1)=cc1;//Torna o valor inicial de Y igual a cc 1
y(m)=cc2; //Torna o valor final de Y igual a cc 2
x(i)=Mx(i) //torna x iagual ao vetor de pontos do domínio - Mx
y(i)=D*sin(E*x(i))+F*cos(G*x(i))+H*x(i)^3+I*x(i)^2+J*x(i)+L //Igualdade da EDO - f(x) end
ML=MA\y
//DETRMINA A FUNÇÃO Fx SEM AS CONDIÇÕES DE CONTORNO
for i=1:m //Cria um laço cosiderando a função no contorno da fronteira
x(i)=Mx(i) //torna x iagual ao vetor de pontos do domínio - Mx
104
fx(i)=D*sin(E*x(i))+F*cos(G*x(i))+H*x(i)^3+I*x(i)^2+J*x(i)+L //Igualdade da EDO - f(x)
end
//---------CALCULA O RESÍDUO PARA DETERMINAR O "C" ÓTIMO (PF)
RES=sum(abs((LFI*ML-fx))) //Resíduo para cada passo k
MPF(k)=PF //Variável que armazena os valores de C
MSR(k)=RES //Variável que armazena os resíduos
//---------CALCULA O RESÍDUO PARA AS CONDIÇÕES DE CONTORNO
if Escolhacc1==1 & Escolhacc2 ==1 then
RESCC1=sum(abs((FI*ML-cc1))); //armazena os resíduos da primeira cc (cc=u(x)) RESCC2=sum(abs((FI*ML-cc2))); //armazena os resíduos da segunda cc (cc=u(x))
elseif Escolhacc1==1 & Escolhacc2 ==2 then
RESCC1=sum(abs((FI*ML-cc1))); RESCC2=sum(abs((LFI1*ML-cc2)));
elseif Escolhacc1==2 & Escolhacc2 ==2 then RESCC1=sum(abs((LFI1*ML-cc1))); //armazena os resíduos da primeira cc (cc=u'(x))
RESCC2=sum(abs((LFI1*ML-cc2))); //armazena os resíduos da segunda cc (cc=u'(x))
elseif Escolhacc1==2 & Escolhacc2 ==1 then RESCC1=sum(abs((LFI1*ML-cc1))); //armazena os resíduos da primeira cc (cc=u'(x))
RESCC1=sum(abs((FI*ML-cc2)));
end
RESCC=RESCC1+RESCC2; MRC(k)=RESCC //Armazena os residuos das duas condições de contorno
end
//if PF==PFfinal then [MenorResiduo LinhaMSR]=min(MSR) //Localiza o valor do menor resíduo calculado e sua
respectiva linha
PFindicado=MPF(LinhaMSR) //Armazena o valor de C ótimo referente a linha do menor resíduo PF=PFindicado //Parâmetro de forma otimizado
end
if PFindicado<>0 then PF=PFindicado
//----- LAÇO PARA PROCURA DA SOLUÇÃO NUMÉRICA ------
//------------DETERMINA A MATRIZ R
[m n]=size(Mx) //obtem a quantidade de linha (m) e colunas (n) for i=1:m //Laço para linhas
for j=1:m //Laço para as colunas
R(i,j)=([Mx(i,1)-My(j,1)]^2)^0.5; //Determina a matriz R
end
end
//------------DETERMINA A MATRIZ FI
[m n]=size(R) //obtem a quantidade de linha (ri) e colunas (rj)
for i=1:m //Laço para linhas
for j=1:n //Laço para as colunas
105
if Metodo == 2 then
FI(i,j)=(R(i,j)^2+PF^2)^(-0.5);
else FI(i,j)=(R(i,j)^2+PF^2)^0.5; //Determina a matriz FI
end//end if
end //end for j
end// end for i
//------------CALCULA A MATRIZ L(FI)
[m n]=size(FI);
//Determina L(FI) para 2ª derivada for i=1:m //Laço para linhas
for j=1:n //Laço para as colunas
if Metodo == 2 then
RFI2(i,j)=(R(i,j)^2+PF^2)^(-2.5); LFI2(i,j)=RFI2(i,j)*((-(R(i,j)^2+PF^2))+(3*R(i,j)^2));//Matriz L(FI)
else
RFI2(i,j)=(R(i,j)^2+PF^2)^1.5; LFI2(i,j)=PF^2/RFI2(i,j); //Matriz L(FI)
end //end if
end//end for j
end//end for i
//Determina L(FI) para 1ª derivada
for i=1:m //Laço para linhas
for j=1:n //Laço para as colunas if Metodo == 2 then
LFI1(i,j)=-R(i,j)/(R(i,j)^2+PF^2)^(1.5);
else
LFI1(i,j)=R(i,j)/(R(i,j)^2+PF^2)^0.5; end//end if
end//end for j
end//end for i
//-------------MONTA A MATRIZ L(FI)
LFI=(A*LFI2)+(B*LFI1)+(C*FI); //Determina a matriz LFI
//-------------DETERMINA A MATRIZ 'A' E SUA INVERSA
if Escolhacc1 == 2 then
A1=LFI1(1,1:n);
else A1=FI(1,1:n);//recorta a 1ª linha da matriz FI (cc1)
end
//A1=FI(1,1:n); //recorta a 1ª linha da matriz FI (primeira condição cont.)
A2=LFI(2:m-1,1:n); //recorta na matriz LFI (pontos interior do domínio)
if Escolhacc2 == 2 then//avalia se deve usar u ou du para cc2
A3=LFI1(m,1:n); //recorta a LFI1 (cc 2 "A derivada de FI") else
A3=FI(m,1:n);//recorta a última linha da FI (cc 2 u(x2)=u2.)
end
106
MA=[A1; A2; A3]; //Concatena linhas para formar a matriz A
//if det(MA)<>0 then
// MAinv=inv(MA); //else
// disp('A matriz não possui inversa') //Avisa se a matriz não for invertível
//abort
//end
//-----------DETERMINA O VETOR Y
[m n]=size(Mx);
y(1)=cc1;//Torna o valor inicial de Y igual a cc 1 y(m)=cc2; //Torna o valor final de Y igual a cc 2
for i=2:m-1 //Cria um laço desconsiderando as cc na fronteira
x(i)=Mx(i); //torna x iagual ao vetor de pontos do domínio - Mx
//Igualdade da EDO - f(x) y(i)=D*sin(E*x(i))+F*cos(G*x(i))+H*x(i)^3+I*x(i)^2+J*x(i)+L;
end
//---------DETERMINA A MATRIZ LAMBDA - ML
ML=MA\y; //produto da matriz A inversa pelo vetor y
//DETERMINA A MATRIZ DE APROXIMAÇÃO - Map = FI*ML Map=FI*ML; //matriz com os pontos de aproximação
end
//APRESENTA OS RESULTADOS EM GRÁFICO
clf(); //Limpa a tela gráfica
xgrid(2) //Coloca uma grade no gráfico
plot(Mx,Map,'o-b'); //Plota um gráfico com os pontos aproximado
//SOLUÇÃO ANALÍTICA DOS ALUNOS DA TURMA INSERIDOS NAS LINHAS SEGUINTES
x=Mx
z =-exp(-3*x)/[3*(exp(3)-exp(-3))]+exp(3*x)/[3*(exp(3)-exp(-3))]-x/3;
plot(x,z,'ro'); //PLOTA A SOLUÇÃO ANALÍTICA
xlabel(["Pontos do Domínio - x"]);
ylabel("Valor Aproximado ; Valor Exato - u(x)");
title("Gráfico da EDO"); legend(["Solução Numérica - MESHLESS";"Solução Exata"]);
disp(Map); //Mostra na tela o pontos de aproximação
if MSR<>0 then scf()
plot(MPF,MSR,'ro-'); //Plota residuo
xgrid(2) //Coloca uma grade no gráfico xlabel(["Parâmetro de Forma - C"]);
ylabel("Resíduo");
title("Gráfico de Resíduo");
txt=string(PFindicado);
107
leg=legend("Parâmetro de Forma Ótimo no Domínio: " +[txt],1);
end
//OBTEM O VELOR MINIMO DE MSC
[MenorResiduoCC LinhaMRC]=min(MRC) //Localiza o valor do menor resíduo calculado e sua
respectiva linha
PFCC=MPF(LinhaMRC) //Armazena o valor de C ótimo referente a linha do menor resíduo txt1=string(PFCC);
if MRC<>[] then
scf() plot(MPF,MRC,'o-b'); //Plota residuo das cc
xgrid(2) //Coloca uma grade no gráfico
xlabel(["Parâmetro de Forma - C"]);
ylabel("Resíduo das Condições de Contorno"); title("Gráfico de Resíduo nas Condições Contorno e no Domínio");
//legend(["Parâmetro de Forma no Contorno: "]);
leg=legend("Parâmetro de Forma Ótimo no Contorno(C): " +[txt1],1); end
//FIM
108
APÊNDICE E- ROTINA DO PROBLEMA DA EQUAÇÃO DE DIFUSÃO-PURA. // MESTRANDO: FRANCISCO REIS
clear;clf;
//Descrição do problema
//Equação diferencial: ∂²u/∂t =D.∂²u/∂x² //Condição inicial: x(x,0)=f(x)=2x(1-x)
//Condição de contorno: u(0,t)=u(L,t)=0
//função para a solução analítica function y=u(x, t)
soma=0;
for n=0:1000
a=2*n+1; b=16/(a*%pi)^3;
c=exp(-2*a^2*%pi^2*t);
d=sin(a*%pi*x);
A=b*c*d; soma=soma+A;
end
y=soma;
endfunction //Função para criar os vetores no espaço e no tempo
function [x, t]=vetores(x0, xf, t0, tf, N)
x=linspace(x0,xf,N); t=linspace(t0,tf,N);
endfunction
//Função meshless/////////////////////////////////////////////////////// function [r, FI, DFI, D2FI, LFI, CI, G, A, LAMBDA]...
=meshless(x,N,c)
//matrizes do meshless
for i=1:N for j=1:N
r(i,j)=x(i)-x(j);
FI(i,j)=sqrt(r(i,j)^2+c^2);
DFI(i,j)=r(i,j)/FI(i,j); D2FI(i,j)=c^2/(FI(i,j)^3)
end
end //matriz LFI
LFI_1=FI;
LFI_2=(-D*dt)*D2FI;
LFI=LFI_1+LFI_2; //vetor condição inicial
for i=1:N
CI(i)=2*x(i)*(1-x(i));
end //concatenação da condição inicial e de contorno c(0,t)=c(L,t)=0
G=[0;CI;0];
//Matriz de colocação A=[FI(1,:);LFI;FI(N,:)];
//Cálculo do lambda
109
LAMBDA=A\G;
endfunction
//Dados do problema x0=0;xf=1;
t0=0;tf=0.20;
D=2;//Coeficiente de dispersão
NA=200;//Número de pontos da solução analítica; NN=40; //Número de pontos da solução numérica;
//Solução analítica
t0=0*tf;
t1=0.25*tf; t2=0.50*tf
[x,t]=vetores(x0,xf,t0,tf,NA);
ANALITICA_t0=u(x,t0);
ANALITICA_t1=u(x,t1); ANALITICA_t2=u(x,t2);
plot(x,ANALITICA_t0,x,ANALITICA_t1,x,ANALITICA_t2); xlabel("$x(Km)$","fontsize", 4);
ylabel("$c(x,t)$","fontsize", 4);
hl=legend(['t=0s';'t=0.05s';'t=0.10s']);
hl=legend(['Num.t=0s';'Num.t=0.05s';'Num.t=0.10s']); //solução numérica
[x,t]=vetores(x0,xf,t0,tf,NN);
dx=x(2)-x(1);
dt=(dx^2)/6*D;
for j=1:100
c(j)=j/100;
[r,FI,DFI,D2FI,LFI,CI,G,A,LAMBDA_0]=meshless(x,NN,c(j)); VR=A*LAMBDA_0-G;
RES=sum(VR^2);
R_AC(j)=RES; end
[R_min, pos_c_otimo]=min(R_AC);
printf("C ótimo = %g", c(pos_c_otimo)); [r,FI,DFI,D2FI,LFI,CI,G,A,LAMBDA_0]=meshless(x,NN,c(pos_c_otimo));
NUMERICA_t0=FI*LAMBDA_0;
tempo=0; k=0;
while tempo<=tf
k=k+1;
tempo=tempo+dt; sol=FI*LAMBDA_0;
B=[0;sol;0];
LAMBDA=(A\B); sol=FI*LAMBDA;
for i=1:NN
NUMERICA(k,i)=sol(i);
end
110
LAMBDA_0=LAMBDA;
end
NUMERICA_t1=NUMERICA(int(t1/dt),:);
NUMERICA_t2=NUMERICA(int(t2/dt),:);
plot(x,NUMERICA_t0,"*",x,NUMERICA_t1,"*",x,NUMERICA_t2,"*");
111
APÊNDICE F- ROTINA DO PROBLEMA DA EQUAÇÃO DE DIFUSÃO-
ADVECÇÃO PARA OS TEMPOS DE 12 min, 24 min e 36 min.
clear;
stacksize("max");
//Dados do problema*************************************************** z=1;L=60;V=2.4;D=0.15;m=5000;
Area=z*L;mp=m/Area;
//Dominio*************************************************************
x0=0; xl=1;
xf=6;
//********************************************************************
t1=0.2;t2=2*t1;t3=3*t1; dx=0.09;
dt=0.05*dx/V;
Courant=V*dt/dx; N=round((xf-x0)/dx)+1;
N1=round((xl-x0)/dx)+1;
xi=[linspace(x0,xl,N1) linspace(xl+dx,xf,N-N1)];
xj=[linspace(x0,xl,N1) linspace(xl+dx,xf,N-N1)]; [oi N]=size(xi);
//Função Meshless*****************************************************
//function [M,R]=meshless(c,xi,xj,N,V,D,xl,t1,t2,t3)
function [M, R_media, LAMBDA, DLT]=meshless(c, xi, xj, N, V, D, xl, t1, t2, t3) //inicializar as matrizes vazias**************************************
r=[];FI=[];DFI=[];D2FI=[];LFI=[];F1=[];A=[];B=[];M=[];
//matrizes do meshless************************************************
for i=1:N for j=1:N
r(i,j)=xi(i)-xj(j);
end end
FI=sqrt(r.^2+c^2);
DFI=r.*(FI).^(-1);
D2FI=c^2*(FI.^(-3)); //matriz LFI**********************************************************
//D=D+(V^2*dt)/2;//Dispersão numérica (livro do chapra)
LFI=FI+dt*(V*DFI-D*D2FI);
//vetor condições iniciais******************************************** //A concentração é nula em todos os pontos exceto no lançamento x=xl**
for i=1:N
x=xi(i)
if x==xl then F1(i)=mp/(dx*1000)
else
F1(i)=0; end
end
//********************************************************************
//vetor condições de contorno***************************************** B(1)=0;//condição de contorno
112
B(N)=0;//condição de contorno
//********************************************************************
LAMBDA_ANT=FI\F1 ///????????????????????? tempo=0;
k=0;
A=LFI;
A(1,:)=FI(1,:);//matriz de colocação A(N,:)=FI(N,:);//matriz de colocação
//laço do tempo*******************************************************
while (tempo<t3)
tempo=tempo+dt; k=k+1;
FIL=FI*LAMBDA_ANT
B(2:N-1)=FIL(2:N-1);
LAMBDA=A\B;// LAMDA=B|F1 DLT=(LAMBDA-LAMBDA_ANT)/dt;
//R=V*DFI*LAMBDA-D*D2FI*LAMBDA;// //vetor resíduo
R=FI*DLT+V*DFI*LAMBDA-D*D2FI*LAMBDA; //R_ANT=FI*DLT+V*DFI*LAMBDA_ANT-D2FI*LAMBDA_ANT;
//R=R+R_ANT;
// LAMBDA_ANT=LAMBDA
// SOL=FI*LAMBDA if (k==1) then
R_AC=(abs(R))
else
R_AC=R_AC+(abs(R)) end
SOL=FI*LAMBDA;
LAMBDA_ANT=LAMBDA;
for i=1:N M(k,i)=SOL(i)
end
end R_media=sum(R_AC/k)
//********************************************************************
endfunction
//******************************************************************** //função para a solução analítica*************************************
function y=C(t, x)
C1=mp/(sqrt(4*%pi*D*t)*1000);
C2=(x-xl-V*t)^2; C3=4*D*t;
y=C1*exp(-C2/C3);
endfunction
//******************************************************************** c_0=0.01;
c_f=1;
d_c=0.01; N_INT = int((c_f - c_0)/d_c)+1
Ds=xf-xl;
NS=N;
113
dc=Ds/(NS-1);
//c0=0.01; //c=linspace(0.001,1,100);
//c=[];
winH=waitbar("Aguarde...")
waitbar(0,winH) nint=20
for i=1:nint
alfaS=2+(i-1)*d_c;
c(i)=alfaS*dc//rand(); //c(i)=c_0 + (i-1)*d_c;
[M,R_media,LAMBDA,DLT]=meshless(c(i),xi,xj,N,V,D,xl,t1,t2,t3);
R_AC(i)=R_media;
waitbar(i/nint,winH) end
close(winH);
//for i=1:N_INT
// //alfaS=2.0+(i-1)*d_c;
// //c(i)=alfaS*dc//rand();//c_0 + (i-1)*d_c;
// c(i)=c_0 + (i-1)*d_c; // [M,R,LAMBDA]=meshless(c(i),xi,xj,N,V,D,xl,t1,t2,t3)
// R_AC(i)=sum(abs(R));
//end
subplot(211);
//clf(1);
plot(c,R_AC);
xlabel("$c$","fontsize", 4); ylabel("$R$","fontsize", 4);
[R_min,pos_c_otimo]=min(R_AC);
printf("c ótimo = %g\n",c(pos_c_otimo));
//Passo 4: Plota a solução analítica e a numérica para c=c_otimo
//********************************************************************
//Chamada do meshlees com a otimização do c*************************** [M,R,LAMBDA,DLT]=meshless(c(pos_c_otimo),xi,xj,N,V,D,xl,t1,t2,t3);
//*******************************************************************
//Plotagem dos gráficos***********************************************
subplot(212); //clf(2)
XAN=[linspace(x0,xl,50) linspace(xl,xf,300)]
plot(XAN,C(t1,XAN),"-b",XAN,C(t2,XAN),"-g",XAN,C(t3,XAN),"-
r",xi,M(t1/dt,:),".b",xi,M(t2/dt,:),".g",xi,M(t3/dt,:),".r"); //plot(xi,C(t1),xi,C(t2),xi,C(t3));
xlabel("$x\,(\text{km})$","fontsize", 4);
ylabel("$c\,(\text{mg}/\text{L})$","fontsize", 4); hl=legend(['t=12min';'t=24min';'t=36min';'Num.t=12min';'Num.t=24min';'Num.t=36min']);
//Problema//Um rio tem os seguintes dados:
//Profundidade z=1m;
//Largura L=60m;
114
//Velocidade V=2400 m/h=2.4 km/h;
//Dispersão D=150000 m2/h=0.15km^2/h;
//Um derramamento de uma substância conservativa no valor de 5kg=5000g //ocorre em um rio a 0.5km abaixo de uma cidade.
//problema formulado:
//Equação diferencial: dc/dt+V*dc/dx=D*d^2c/dx^2;
//solução analítica: //c(x,t)=(mp/sqrt(4*%pi*D*t))*exp(-(x-xl-V*t)/(4*D*t));
//Condições iniciais:
//c(x,0) é nula em todos os pontos, exceto sobre a seção reta definida
//por x=xl, onde foi introduzido o poluente. Esta afirmação pode ser //descrita por meio da função delta de dirac
//Condições de contorno:
//A concentração é nula nos seguintes pontos:
//c(x,t)=c(+inf,t)=c(-inf,t)=0 //********************************************************************
//Otimizaçao do c*****************************************************
//r_menor=-10^25 //c0=(1/1000);
//for (j=0:1:10)
// //c=(2*rand())/sqrt(N);
// c=c0*(j+1); // [M,R]=meshless(c,xi,xj,N,V,D,xl,t1,t2,t3)
// r_atual=sum(abs(R))
// if r_menor<=r_atual then
// r_menor=r_atual // c_menor=c
// end
//end
//******************************************************************** //c0=0.001;
//[M,R]=meshless(c0,xi,xj,N,V,D,xl,t1,t2,t3) ;
//r_atual=sum(abs(R)) //r_menor=+10^250;
//k=0;
//c_menor=2/sqrt(N);
//c=c_menor; //while k<=10 do
// k=k+1;
// //c=c0*(1+k);
// // [M,R]=meshless(c,xi,xj,N,V,D,xl,t1,t2,t3)
// r_atual=sum(abs(R));
// if r_atual<=r_menor then
// r_menor=r_atual // c_menor=c;
// end
// // printf("k=%g\n",k)
// if k==50 then
// break
// end/
115
APÊNDICE G - ROTINA DA APROXIMAÇÃO DA SOLUÇÃO NUMÉRICA COM A
SOLUÇÃO ANALÍTICA DO PROBLEMA DA EQUAÇÃO DE DIFUSÃO-ADVECÇÃO CONSIDERANDO O TEMPO DE 15 min.
clear; //Dados do problema******************************************************
z=1;//profundidade (m)
L=60;//largura(m)
Area=z*L;//Area(m^2) V=2.4;//velocidade (km/h)
D=0.15;//Dispersão em (m^2/h)
m=5*10^3;//massa do poluente (g)
mp=m/Area;//massa por unidade de área. (g/m2) //***********************************************************************
//função para a solução analítica****************************************
function [y]=C(x, t);
C1=(0.001*mp)/(sqrt(4*%pi*D*t)); C2=(x-xL-V*t)^2;
C3=4*D*t;
y=C1*exp(-C2/C3); //Obs.: Multiplica-se por 0.001 para converter as unidades!
endfunction
//***********************************************************************
//Função para criar os vetores x e t************************************* function [x, t, dx, dt, xL]=vetores(x0, xf, t0, tf, N)
x=linspace(x0,xf,N); t=linspace(t0,tf,N);
dx=x(2)-x(1);
dt=(Courant*dx)/V; [a b]=size(x);
p=int(b/5);
xL=x(p);
endfunction //***********************************************************************
//Função meshless
function [r, FI, DFI, D2FI, LFI, CI, G, A, LAMBDA]... =meshless(x,N,c)
//matrizes do meshless **************************************************
for i=1:N
for j=1:N r(i,j)=x(i)-x(j);
FI(i,j)=sqrt(r(i,j)^2+c^2);
DFI(i,j)=r(i,j)/FI(i,j);
D2FI(i,j)=c^2/(FI(i,j)^3) end
end
//matriz LFI*************************************************************
LFI_1=FI; LFI_2=V*dt*DFI;
LFI_3=-D*dt*D2FI;
LFI=LFI_1+LFI_2+LFI_3; //vetor condição inicial*************************************************
//A concentração é nula em todos os pontos exceto no lançamento x=xL
116
for i=1:N
if x(i)==xL then
CI(i)=(0.001*mp)/dx; else
CI(i)=0;
end
end //***********************************************************************
//concatenação da condição inicial e de contorno c(x0,t)=c(xf,t)=0*******
G=[0;CI;0];
//Matriz de colocação**************************************************** A=[FI(1,:);LFI;FI(N,:)];
//Cálculo do lambda
LAMBDA=A\G;
endfunction //fim da função**********************************************************
//Dados de entrada dos espaços e dos tempos******************************
x0=0;//posição inicial em km xf=5;//posição final em km
t0=0;//tempo inicial em h
tf=1;//tempo final em h.
Courant=0.01;//número de courant t1=0.25*tf;t2=2*t1;t3=3*t1;//tempos para plotagem dos gráficos
N_analitica=300;//número de pontos da solução analítica.
N_numerica=100;//número de pontos para a solução numérica.
c=0.15;//rand();//parâmetro de forma do meshless //Solução analítica******************************************************
[x,t,dx,dt,xL]=vetores(x0,xf,t0,tf,N_analitica);
ANALITICA_t1=C(x,t1);ANALITICA_t2=C(x,t2);ANALITICA_t3=C(x,t3);
plot(x,ANALITICA_t1); //plot(x,ANALITICA_t1,x,ANALITICA_t2,x,ANALITICA_t3);
xlabel("$x\,(\text{km})$","fontsize", 4);
ylabel("$c\,(\text{mg}/\text{L})$","fontsize", 4); hl=legend(['t=15 min';'t=30 min';'t=45 min']);
//***********************************************************************
//Solução numérica*******************************************************
[x,t,dx,dt,xL]=vetores(x0,xf,t0,tf,N_numerica);
[r,FI,DFI,D2FI,LFI,CI,G,A,LAMBDA_0]...
=meshless(x,N_numerica,c)
//Laço para o tempo***************************************************** tempo=0;
k=0;
while tempo<=t1
k=k+1; tempo=tempo+dt;
sol=FI*LAMBDA_0;
B=[0;sol;0]; LAMBDA=(A\B);
VR=A*LAMBDA-B;
R(k)=sum(VR)^2;
sol=FI*LAMBDA;
117
for i=1:N_numerica
NUMERICA(k,i)=sol(i);
end LAMBDA_0=LAMBDA
end
TN1=NUMERICA(int(t1/dt),:);
//TN2=NUMERICA(int(t2/dt),:);
//TN3=NUMERICA(int(t3/dt),:);
plot(x,TN1,"."); //plot(x,TN1,x,TN2,x,TN3);
