TT009 Matemática Aplicada I
Curso de Engenharia Ambiental
Departamento de Engenharia Ambiental, UFPR
Trabalho Computacional. Data de entrega: 01 mai 2012
Prof. Nelson Luís Dias
1
Instruções para a entrega do trabalho
Atenção: Você entregará apenas 2 arquivos digitais, que você deverá enviar para
[email protected],
com o assunto: TT009-tc1-xxxx . Os nomes dos arquivos serão
1. rksub-xxxx.py
2. rksub-xxxx.pdf
Nos nomes acima, xxxx é o seu código individual de identificação na disciplina, disponível na tabela 1. O seu programa deve estar escrito em Python ≥ 2.7 e < 3.x. O
arquivo rksub-xxxx.py deve utilizar codificação ISO-8859-1; o programa deve gerar
obrigatoriamente apenas um arquivo de saída, cujo nome será
rksub-xxxx.dat
O arquivo rksub-xxxx.pdf deve conter suas respostas ao enunciado mais abaixo,
figuras, equações utilizadas, referências bibliográficas, etc.. Não há regras muito fixas,
mas você deve obedecer ao seguinte:
• 2,5 cm de margem (direita, esquerda, acima e abaixo)
• espaçamento simples
• Tipo romano com serifa (por exemplo, Times Roman ou Times New Roman; não
use tipos sem serifa, tais como Arial ou Calibri)
O arquivo rksub-xxxx.dat deve ser um arquivo texto na codificação ISO-8859-1;
ele deve ser organizado em 2 colunas de 10 caracteres cada, impressas com o formato
%10.6f, sem brancos adicionais entre as colunas. A primeira coluna deve conter os valores de ξ , a partir de 0, em incrementos ∆ξ que você definirá na sua solução; a segunda
coluna deve conter os valores de φ (ξ ) (vide seção 4).
2
Correção do trabalho
• Se dois ou mais trabalhos tiverem elementos comuns que configurem cópia, a nota
destes trabalhos é zero.
• Se o seu programa não rodar, sua nota é zero.
• Se o seu programa não gerar os arquivos de saída com os nomes e o formato especificados, sua nota é zero.
• Eu vou comparar seu arquivo rksub-xxxx.dat contendo a sua solução numérica
com a solução numérica correta.
1
Quando seus arquivos passarem nos testes acima, então eu vou verificar seu código
em rksub-xxxx.py e ler seu trabalho em rksub-xxxx.pdf, para então dar a nota final.
A qualidade da apresentação, a correção do Português, e a obediência às regras
usuais de redação acadêmica (por exemplo, uso correto de citações e referências
bibliográficas) serão fortemente levados em consideração na sua nota.
2
Tabela 1: Códigos de identificação individuais
Nome
xxxx
ALINE MEDEIROS FERREIRA DE ARAÚJO
ANDRE LUIZ FREIRE DE CARVALHO KATO
ANGELICA MARTINI
BARBARA DE TOLEDO MONTANDON DUMONT
BARBARA WOLFF ZWOLINSKI
DIANDRA CHRISTINE VICENTE DE LIMA
DIEGO TAKESHI MIYASAKA
EINARA ZAHN
EVELIN DE LARA PALLU
FABIO BORTOLOTO VALEBONA
FRANCIELE STELMACHUK
GABRIEL CIRINO SAMY
GABRIELA GHOBAD EL SAYED
GABRIELA GUERRIZE CONTE
GABRIELE STURM
GISMAR FERNANDES RIBEIRO
GUILHERME TANFERRI
HENRIQUE GUARNERI
IRMA CATALINA SALAZAR BAY
ISABELA CHEMIN
ISADORA CHIAMULERA
JESSYCA PETRY DALAZEN
JOAO PEDRO BAZZO VIEIRA
JULIA NAGAFUTI DOS SANTOS
LEON FERNANDO MIECOANSKI
LUÍS FELIPE PANKIVEVICZ
MARCELA ODORIZZI
MARCELO LUIZ NORILLER
MARIANA ANGELITA MOISÉS
MATHEUS HENRIQUE TAVARES
MELINA NAGATA BELTRANE
MURIEL EDYTH LUMSDEN SZYMANSKI PATRICIO
RENA DOMINIKA SALZMANN
THIAGO MORIGGI
VICTOR KYOCHI BERNARDES
arau
kato
mart
dumo
swol
lima
miya
zahn
pall
vale
stel
samy
saye
cont
stur
ribe
tanf
guar
sbay
chem
chia
dala
viei
sant
miec
pank
odor
nori
mois
tava
belt
patr
salz
mori
bern
3
3
Exemplo de um trabalho nota 100
Nesta seção eu dou um exemplo de um trabalho “nota 100”. É claro que o exemplo é mais
simples do que o deste trabalho “real”, mas ele deve lhe dar uma boa idéia de como fazer
as coisas.
Enunciado:
Utilizando o método de Runge-Kutta (!), calcule a função gama
Γ(x) ≡
Z ∞
0
t x−1 e−t dt
numericamente para o intervalo x ∈ [1, 5].
Solução:
A variável independente x entra com um parâmetro, não como um dos limites de
integração. Isto torna as coisas um pouco mais complicadas. Para obter a solução, vamos
primeiro definir a função gama incompleta inferior
(http://en.wikipedia.org/wiki/Incomplete_gamma_function
e Abramowitz e Stegun (1972)):
γ(x, y) ≡
Z y
0
t x−1 e−t dt.
(1)
Toda integração numérica pode ser convertida na solução numérica de uma EDO
(Press et al., 1992, p. 129). Por exemplo, no caso acima, a integração numérica é equivalente a
dγ(x,t)
= t x−1 e−t ;
γ(x, 0) = 0.
(2)
dt
Se resolvermos a equação diferencial acima, então teremos, simplesmente,
Γ(x) = γ(x, ∞).
(3)
Numericamente, nós vamos substituir ∞ por um valor finito porém suficientemente grande.
Para termos uma boa idéia do que signfica um valor “suficientemente grande”, vamos plotar o integrando t x−1 e−t para os valores x = 1, x = 2, x = 3, x = 4 e x = 5. Isto é mostrado
na figura 1.
Note que o limite prático de integração varia com x. No entanto, por simplicidade,
podemos adotar t = 20 para todos os casos (com um pouco de exagero). Em outras
palavras, nós esperamos que (numericamente) Γ(x) ≈ γ(x, 20) para 1 ≤ x ≤ 5.
Um outro ponto importante é que é relativamente simples provar que
Γ(x) = (x − 1)!
(4)
para x ∈ N. Portanto, é imediato que Γ(1) = 1, Γ(2) = 1, Γ(3) = 2, Γ(4) = 6 e Γ(5) = 24.
Isto significa que é possível testar a solução numérica para estes valores particulares de x.
Nosso primeiro esforço, portanto, é escrever um protótipo para testar alguns ou todos
entre estes valores particulares. Lembremo-nos de que isto é feito simplesmente integrando a equação diferencial (2) até 20(?), com um passo que precisa ser verificado por
tentativa e erro.
4
Listagem 1: Primeiro protótipo para o cálculo de Γ(x)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
# !/ usr / bin / python
# -* - coding : iso -8859 -1 -* # ---------------------------------------------------------# intgam0 : resolve a integral
#
# \ int_0 ^ infty t ^{ x -1} e ^{ - t }\ , dt
#
# usando o método de Runge - Kutta de ordem 4
#
# 2012 -02 -14 T20 :45:29
# 2012 -02 -14 T20 :45:40
#
# uso : ./ intgam0 . py
# ---------------------------------------------------------from __future__ import prin t_func tion
from __future__ import division
from math import exp
h = 1.0 e -4
# passo de integração
to = 0.0
# t inicial
go = 0.0
# gamma inicial
n = int (30.0/ h )
# número de passos até o " infinito "
# ---------------------------------------------------------# define a função a ser integrada
# ---------------------------------------------------------x = 5.0
def ff (t , g ):
return t **( x -1.0)* exp ( - t )
def rk4 (x ,y ,h , ff ):
'''
rk4 implementa um passo do método de Runge - Kutta de ordem 4
'''
k1 = h * ff (x , y )
k2 = h * ff ( x + h /2 , y + k1 /2)
k3 = h * ff ( x + h /2 , y + k2 /2)
k4 = h * ff ( x +h , y + k3 )
yn = y + k1 /6.0 + k2 /3.0 + k3 /3.0 + k4 /6.0
return yn
for i in range (0 , n ):
# loop da solução numérica
tn = ( i +1)* h
# novo t
gn = rk4 ( to , go ,h , ff )
# novo gamma (x , t )
( to , go ) = ( tn , gn )
# o antigo passa a ser o novo
# ---------------------------------------------------------# imprime o resultado final
# ---------------------------------------------------------print ( ' %10.6 f ␣ %10.6 f ' % ( tn , gn ) )
5
5
x=1
x=2
x=3
x=4
x=5
dγ(x,t)/dt
4
3
2
1
0
0
2
4
6
8
10
12
14
16
18
20
t
Figura 1: O integrando γ(x,t) para x = 1, . . . , 5.
Nosso primeiro protótipo, muito simples, e denominado intgam0.py, é mostrado na
listagem 1. Note que nós reaproveitamos a rotina rk4 que foi apresentada em sala de aula.
Nós observamos, por tentativa e erro, que h = 10−4 é um passo suficientemente “preciso” para o método de Runge-Kutta, mas que o “infinito” é melhor representado numericamente por 30, e não por 20. Observe também na listagem 1 que o valor de x é definido
“na força bruta”, de maneira um pouco deselegante, por meio de uma variável global dentro do próprio programa. Para testar x = 1, . . . , 5, nós simplesmente editamos o programa
e modificamos o valor de x com um editor de texto.
O programa “definitivo” chama-se intgam1.py. Ele é mostrado na listagem 2. A
principal diferença em relação a intgam0 consiste em um loop externo, implementado
por meio de um while, em que o valor de x, que antes era fixo, é incrementado de 0,01
em 0,01. Além disto, os valores para cada x são impressos em um arquivo de saída.
O resultado de intgam1.py é mostrado na figura 2. Na figura, note que a concordância é (visualmente) perfeita. Algumas observações importantes são as seguintes:
• As figuras deste trabalho foram geradas com o programa de plotagem Gnuplot
(www.gnuplot.info). Como Gnuplot possui uma função gamma(x) “embutida”,
ela foi utilizada para comparação com os pontos calculados por intgam1.py na
figura 2.
• Apenas 1 a cada 5 pontos gerados está plotado na figura 2 para não sobrecarregar a
figura.
Conclusões
Neste trabalho, nós mostramos que é possível utilizar um método de solução numérica
de equações diferenciais ordinárias (o método de Runge-Kutta de 4a ordem) para integrar
numericamente uma função. O teste foi realizado com a função Γ(x). Devido ao fato
de que x neste problema é um parâmetro dentro do integrando, e não um dos limites
6
Listagem 2: Programa definitivo para o cálculo de Γ(x)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
# !/ usr / bin / python
# -* - coding : iso -8859 -1 -* # ---------------------------------------------------------# intgam0 : resolve a integral
#
# \ int_0 ^ infty t ^{ x -1} e ^{ - t }\ , dt
#
# usando o método de Runge - Kutta de ordem 4
#
# para x =1 ,1+ dx ,... ,5 com dx = 0.01
#
# 2012 -02 -14 T21 :16:09
# 2012 -02 -14 T21 :16:14
#
# uso : ./ intgam1 . py
# ---------------------------------------------------------from __future__ import prin t_func tion
from __future__ import division
from math import exp
h = 1.0 e -4
# passo de integração de rk4
dx = 0.01
# passo de variação de x
n = int (30.0/ h )
# número de passos até o " infinito "
# ---------------------------------------------------------# define a função a ser integrada
# ---------------------------------------------------------def ff (t , g ):
return t **( x -1.0)* exp ( - t )
def rk4 (x ,y ,h , ff ):
'''
rk4 implementa um passo do método de Runge - Kutta de ordem 4
'''
k1 = h * ff (x , y )
k2 = h * ff ( x + h /2 , y + k1 /2)
k3 = h * ff ( x + h /2 , y + k2 /2)
k4 = h * ff ( x +h , y + k3 )
yn = y + k1 /6.0 + k2 /3.0 + k3 /3.0 + k4 /6.0
return yn
fou = open ( ' intgam1 . dat ' , ' wt ')
# abre o arquivo de saída
x = 1.0
while ( x <= 5.0) :
to = 0.0
# t inicial
go = 0.0
# gamma inicial
for i in range (0 , n ):
# loop da solução numérica
tn = ( i +1)* h
# novo t
gn = rk4 ( to , go ,h , ff )
# novo gamma (x , t )
( to , go ) = ( tn , gn )
# o antigo passa a ser o novo
# ---------------------------------------------------------# imprime o resultado no arquivo de saída
# ---------------------------------------------------------print (x , go )
# imprime x na tela
fou . write ( ' %10.6 f ␣ %10.6 f \ n ' % (x , gn ) )
x += dx
fou . close ()
7
25
Γ(x) — intgama1.py
Γ(x) — Gnuplot
20
Γ(x)
15
10
5
0
1
2
3
4
5
x
Figura 2: Cálculo de Γ(x) com o método de Runge-Kutta: saída de intgam1 comparada
com a função Γ(x) do programa de plotagem (Gnuplot)
de integração, o cálculo é mais complicado do que neste último caso, e exige que nós
reconheçamos, antes de mais nada, a necessidade de utilizar uma função gama incompleta
inferior.
8
4
Trabalho computacional
Considere a equação diferencial ordinária
d
dφ
dφ
= 0.
φ
+ 2ξ
dξ
dξ
dξ
(5)
com condições de contorno
φ (0) = 0,
φ (∞) = 1.
Resolva (5) utilizando o método de Runge-Kutta de 4a ordem, conforme esboçado
no material didático distribuído, na seção 2.8. Estude cuidadosamente a seção antes de
iniciar sua solução numérica. Você deve transformar (5) em um sistema de 2 equações
diferenciais ordinárias de ordem 1, e resolver o passo inicial com o método implícito
descrito na seção 2.8. Somente após este passo você deve utilizar o método de RungeKutta.
Referências
Abramowitz, M. e Stegun, I. A., editores (1972). Handbook of mathematical functions.
Dover Publications, Inc., New York.
Press, W. H., Teukolsky, S. A., Vetterling, W. T., e Flannery, B. P. (1992). Numerical
Recipes in C. Cambridge University Press, Cambridge, UK. 1020 pp.
9