Matemtica Aplicada EngenhariaNelson Lus Dias Departamento de Engenharia Ambiental e Lemma Laboratrio de Estudos em Monitoramento e Modelagem Ambiental Universidade Federal do Paran 30 de setembro de 2011

Sumrio

2

Sumrio1 Ferramentas computacionais 1.1 Antes de comear a trabalhar, . . . . . . . . . . . . . . . 1.2 Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Strings e Inteiros . . . . . . . . . . . . . . . . . . 1.2.2 Nmeros de ponto utuante reais e complexos 1.2.3 Obteno de uma curva de permanncia . . . . . 1.2.4 Arquivos texto e arquivos binrios . . . . . . . . . 1.3 Gnuplot . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Ajuste de uma funo . . . . . . . . . . . . . . . 1.4 Python de novo: projeto probabilstico . . . . . . . . . . 1.5 Maxima . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Um 2.1 2.2 2.3 11 11 13 13 15 16 19 22 24 26 30

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

pouco de polinmios, integrais, sries . . . 33 Integrao numrica: motivao . . . . . . . . . . . . . . . . . . . . . 33 A regra do trapzio . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Aproximao de integrais com sries: a funo erro . . . . . . . . . . 42 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 49 50 54 54 61 61 71 83

3 Soluo numrica de equaes diferenciais ordinrias 3.1 Soluo numrica de equaes diferenciais ordinrias . 3.1.1 Soluo numrica; mtodo de Euler . . . . . . . 3.1.2 Um mtodo implcito, com tratamento analtico 3.1.3 Runge-Kutta . . . . . . . . . . . . . . . . . . .

4 Soluo numrica de equaes diferenciais parciais 4.1 Adveco pura: a onda cinemtica . . . . . . . . . . . . . . . . . . . . 4.2 Difuso pura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Difuso em 2 Dimenses: ADI, e equaes elticas . . . . . . . . . . .

A Dados de vazo mdia anual e vazo mxima anual, Rio dos Patos, 19311999 89

3

Sumrio

4

Lista de Tabelas1.1 2.1 Vazes Mximas Anuais (m3 s1 ) no Rio dos Patos, PR, 19311999 . 16 Estimativas numricas de I e seus erros relativos . . . . . . . . . . . 38

5

Sumrio

6

Lista de Figuras1.1 1.2 1.3 2.1 2.2 2.3 2.4 3.1 3.2 3.3 3.4 3.5 3.6 4.1 4.2 4.3 4.4 4.5 Convenes do terminal postscript eps monochrome Times-Roman 18 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 FDA da vazo mxima anual, Rio dos Patos . . . . . . . . . . . . . . 25 Ajuste de uma FDA Weibull aos dados de vazo mxima anual do Rio dos Patos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Integrao numrica de uma funo . . . . . . . . . . . . . . . . . . . 35 A regra do trapzio, com n = 4 e n = 8 trapzios. . . . . . . . . . . . 38 Funo erf(x) calculada por integrao numrica, com trapepsilon e = 1 106 , versus a erf pr-denida em Gnuplot. . . . . . . . . . 44 Funo erf(x) calculada com srie de Taylor, com erf_1, versus a erf pr-denida em Gnuplot. . . . . . . . . . . . . . . . . . . . . . . . . . 45 Soluo da equao (3.1). . . . . . . . . . . . . . . . . . . . . . . . . 50

Comparao da soluo analtica da equao (3.1) com a sada de sucesso.py, para x = 0,01. . . . . . . . . . . . . . . . . . . . . . . 53 Comparao da soluo analtica da equao (3.1) com a sada de sucesso.py, para x = 0,5. . . . . . . . . . . . . . . . . . . . . . . . 53 Comparao da soluo analtica da equao (3.1) com a sada de sucimp.py, para x = 0,5. . . . . . . . . . . . . . . . . . . . . . . . . 55 Comparao da soluo analtica da equao (3.1) com a sada de euler2.py, para x = 0,5. . . . . . . . . . . . . . . . . . . . . . . . . 56 Comparao da soluo analtica da equao (3.1) com a sada de rungek4.py, para x = 0,5. . . . . . . . . . . . . . . . . . . . . . . . 59 Condio inicial da equao 4.3. . . . . . . . . . . . . . . . . . . . . . 61 Soluo numrica produzida por onda1d-ins.py, para t = 250t, 500t e 750t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Soluo numrica produzida por onda1d-lax.py, para t = 500t, 1000t e 1500t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Soluo numrica produzida pelo esquema upwind, para t = 500t, 1000t e 1500t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 Soluo analtica da equao de difuso para t = 0, t = 0,05, t = 0,10 e t = 0,15. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 7

Sumrio 4.6 Soluo numrica com o mtodo explcito (4.35) (crculos) versus a soluo analtica (linha cheia) da equao de difuso para t = 0, t = 0,05, t = 0,10 e t = 0,15. Apenas 1 a cada 5 pontos da grade numrica so mostrados, para facilitar a comparao com a soluo analtica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Soluo numrica com o mtodo implcito (4.39) (crculos) versus a soluo analtica (linha cheia) da equao de difuso para t = 0, t = 0,05, t = 0,10 e t = 0,15. Apenas 1 a cada 5 pontos da grade numrica so mostrados, para facilitar a comparao com a soluo analtica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Soluo numrica com o mtodo de Crank-Nicholson ( (4.45)) (crculos) versus a soluo analtica (linha cheia) da equao de difuso para t = 0, t = 0,05, t = 0,10 e t = 0,15. Apenas 1 a cada 5 pontos da grade numrica so mostrados, para facilitar a comparao com a soluo analtica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Soluo analtica da equao da difuso bidimensional, para t = 0, t = 0, t = 0,1, t = 0,2 e t = 0,3 . . . . . . . . . . . . . . . . . . . . . 4.10 Soluo numrica da equao da difuso bidimensional com o esquema ADI, para t = 0, t = 0, t = 0,1, t = 0,2 e t = 0,3 . . . . . . . . . . .

8

. 77

. 80

. 83 . 87 . 88

Lista de Listagens1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 1.10 1.11 1.12 1.13 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 2.13 2.14 binint.py exemplo de strings, e inteiros. . . . . . . . . . . . . . . floats.py exemplo de uso de float e complex. . . . . . . . . . . patos-medmax.dat vazes mdia e mxima anuais, Rio dos Patos fqiemp.py clculo de uma FDA emprica . . . . . . . . . . . . . . fqiemp.dat FDA emprica da vazo mxima anual no Rio dos Patos. bintext.py Exemplo de arquivo texto e arquivo binrio . . . . . . writearr.py Escreve um arquivo binrio contendo 3 linhas, cada uma das quais com um array de 10 floats. . . . . . . . . . . . readarr.py Escreve um arquivo binrio contendo 3 linhas, cada uma das quais com um array de 10 floats. . . . . . . . . . . . . . . fiqemp.plt programa para plotar a FDA emprica do Rio dos Patos weibull.plt como ajustar uma FDA analtica (Weibull) aos dados da FDA emprica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . interp.py mdulo com funo para calcular uma interpolao linear. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ppbemp.py projeto probabilstico a partir da distribuio emprica mulambdak.max clculo da mdia da distribuio Weibull em em funo de seus parmetros e k . . . . . . . . . . . . . . . . . . . . . achapol.max Polinmio com propriedades denidas . . . . . . . . Sada de achapol.max . . . . . . . . . . . . . . . . . . . . . . . . . . passaquad.max parbola h(x) = ax2 + bx + c passando por (1, f (1)), (3, f (3)) e (5, f (5)). . . . . . . . . . . . . . . . . . . . . . . Sada de passaquad.max . . . . . . . . . . . . . . . . . . . . . . . . . numint.py Integrao numrica, regra do trapzio . . . . . . . . . quadraver1.py Integrao numrica de f (x) com 8 trapzios . . . numint.py Integrao numrica ineciente, com erro absoluto prestabelecido . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . quadraver2.py Integrao numrica ineciente de f (x) com = 0,0001 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . numint.py Integrao numrica eciente, com erro absoluto prestabelecido . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . quadraver3.py Integrao numrica eciente de f (x) com = 0,000001 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vererf.py Clculo da funo erro por integrao numrica . . . . vererf.plt Plotagem da funo erro por integrao numrica versus a erf(x) pr-denida em Gnuplot . . . . . . . . . . . . . . . . . . Clculo de erf(x) com uma srie de Taylor. . . . . . . . . . . . . . . . vererf1.py Clculo da funo erro com srie de Taylor entre 0 e 3. 9 13 15 17 18 19 20 21 22 23 25 28 29 31 34 34 36 37 39 39 40 40 41 41 43 43 46 46

Sumrio 2.15 vererf1.plt Plotagem da funo erro calculada com srie de Taylor versus a erf(x) pr-denida em Gnuplot . . . . . . . . . . . . . . . 2.16 erfs.py Clculo de erf(x) com srie de Taylor, limitado a no mximo 43 termos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 resolve-eqdif Soluo de uma EDO com Maxima . . . . . . . . 3.2 fracasso.py Um programa com o mtodo de Euler que no funciona 3.3 Sada de fracasso.py . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 sucesso.py Um programa com o mtodo de Euler que funciona . 3.5 sucimp.py Mtodo de Euler implcito . . . . . . . . . . . . . . . . 3.6 euler2 Um mtodo explcito de ordem 2 . . . . . . . . . . . . . . 3.7 rungek4 Mtodo de Runge-Kutta, ordem 4 . . . . . . . . . . . . . 4.1 onda1d-ins.py Soluo de uma onda 1D com um mtodo explcito instvel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 surf1d-ins.py Seleciona alguns intervalos de tempo da soluo numrica para plotagem . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 onda1d-lax.py Soluo de uma onda 1D com um mtodo explcito laxtvel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 surf1d-lax.py Seleciona alguns intervalos de tempo da soluo numrica para plotagem . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 difusao1d-ana.py Soluo analtica da equao da difuso . . . . 4.6 divisao1d-ana.py Seleciona alguns instantes de tempo da soluo analtica para visualizao . . . . . . . . . . . . . . . . . . . . . . . . 4.7 difusao1d-exp.py Soluo numrica da equao da difuso: mtodo explcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 divisao1d-exp.py Seleciona alguns instantes de tempo da soluo analtica para visualizao . . . . . . . . . . . . . . . . . . . . . . . . 4.9 alglin.py Exporta uma rotina que resolve um sistema tridiagonal, baseado em Press et al. (1992) . . . . . . . . . . . . . . . . . . . . . . 4.10 difusao1d-imp.py Soluo numrica da equao da difuso: mtodo implcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 difusao1d-ckn.py Soluo numrica da equao da difuso: esquema de Crank-Nicholson. . . . . . . . . . . . . . . . . . . . . . . . patosmedmax.dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

47 48 49 51 51 52 55 57 58 63 64 68 69 73 74 76 77 80 81 84 89

1Ferramentas computacionaisNeste captulo ns fazemos uma introduo muito rpida a 3 ferramentas computacionais. Python uma linguagem de programao, razoavelmente clssica, e da qual ns vamos explorar apenas uma parte chamada de programao procedural. A escolha de Python deve-se ao fato de ela ser uma linguagem fcil de aprender, e de existirem muitas bibliotecas de rotinas em Python que tornam muitas das tarefas de Matemtica Aplicada fceis de implementar em um computador. Gnuplot uma linguagem para gerar grcos. Ela bastante popular, fcil de usar, e muito exvel. Ns vamos gerar todos os nossos grcos bidimensionais (grcos x y) com Gnuplot. Maxima uma linguagem de processamento simblico. Da mesma maneira que ns faremos contas com Python, ns faremos lgebra com Maxima. Os usos que faremos de Maxima estaro longe de explorar todo o seu potencial. Ns vamos apenas calcular algumas integrais, algumas sries de Taylor, e resolver algumas equaes diferenciais ordinrias; entretanto, a rapidez com que faremos isto justica amplamente o seu uso. Infelizmente, a verso de um programa de computador faz diferena, e s vezes faz muita diferena. As verses que eu utilizei para este texto foram: Python: Gnuplot: Maxima: 2.6.5 4.4 5.21.1

Neste momento, existem duas linhagens de Python convivendo: 2.x, e 3.x. O futuro Python 3.x, e eu procurei usar uma sintaxe to prxima de 3.x quanto possvel. Isto facilitado por comandos to tipo from __future__ import , que voc encontrar diversas vezes ao longo do texto. Se voc (j) estiver usando Python 3.x, remova estes comandos do incio de seus arquivos .py.

1.1 Antes de comear a trabalhar,Voc precisar de algumas condies de funcionamento. Eis os requisitos fundamentais: 11

Matemtica Aplicada 1) Saber que sistema operacional voc est usando.

12

2) Saber usar a linha de comando, ou terminal, onde voc datilografa comandos que so em seguida executados. 3) Saber usar um editor de texto, e salvar seus arquivos com codicao ISO-8859-1. 4) Certicar-se de que voc tem Python instalado. 5) Certicar-se de que voc tem Numpy (um mdulo de Python que deve ser instalado parte, e que vamos utilizar seguidamente) instalado. 6) Certicar-se de que voc tem Gnuplot instalado. 7) Certicar-se de que voc tem Maxima instalada. Ateno! Um editor de texto no um processador de texto. Um editor de texto no produz letras de diferentes tamanhos, no cria tabelas, e no insere guras. Um editor de texto reproduz o texto que voc datilografa, em geral com um tipo de largura constante para que as colunas e espaos quem bem claros. Um editor de texto que vem com Windows chama-se notepad, ou bloco de notas nas verses em Portugus; um excelente substituto chama-se notepad2 (http://www.flos-freeware.ch/notepad2.html). Em Linux, editores de texto simples so o gedit, e o kate. Programadores mais experientes costumam preferir o vim, ou o emacs. Estes dois ltimos possuem verses para os 3 sistemas operacionais mais comuns hoje em dia: Windows, Linux e MacOS. Quando voc estiver praticando o uso das ferramentas computacionais descritas neste texto, suas tarefas invariavelmente sero: 1) Criar o arquivo com o programa em Python, Gnuplot, ou Maxima, usando o editor de texto, e salv-lo em codicao ISO-8859-1. Se no for bvio para voc, procure descobrir como voc pode congurar o seu editor para salvar em ISO-8859-1. Se isto for impossvel para voc, trabalhe normalmente, mas no use acentos: em lugar de impossvel, digite impossivel. 2) Ir para a linha de comando. 3) Executar o programa digitando o seu nome (e no clicando!), possivelmente precedido por python, gnuplot, ou maxima. 4) Vericar se o resultado est correto. 5) Se houver erros, voltar para 1), e reiniciar o processo. Neste texto, eu vou partir do princpio de que todas estas condies esto cumpridas por voc, mas no vou detalh-las mais: em geral, sistemas operacionais, editores de texto e ambientes de programao variam com o gosto do fregus: escolha os seus preferidos, e bom trabalho!

13

1.2 Python Listagem 1.1: binint.py exemplo de strings, e inteiros.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

# !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n a = aafro # a uma string # imprime a na tela print ( a ) print ( len ( a )) # imprime o nmero de caracteres de a i = 2**32 - 1 # i o maior int s / sinal que cabe em 32 bits b = unicode ( i ) # converte i na string c orrespon dente b # imprime b na tela print ( b ) print ( --- a : ) # separa a sada do primeiro for # p / cada c de a , imprime seu valor unicode for c in a : print ( ord ( , c , ) = , ord ( c )) print ( --- b : ) # separa a sada do segundo for # p / cada c de b , imprime seu valor unicode for c in b : print ( ord ( , c , ) = , ord ( c )) print ( unichr (227)) # imprime o caractere unicode no 227

Exerccios Propostos 1.1 Voc j devia estar esperando por isto: o que um sistema operacional? Qual o sistema operacional que voc usa? 1.2 Como saber se Python est instalado?

1.2 PythonPython reconhece os seguintes tipos bsicos de variveis: strings, nmeros inteiros, nmeros de ponto utuante, e nmeros complexos. 1.2.1 Strings e Inteiros Strings, ou cadeias de caracteres, so criaturas do tipo abacaxi , e nmeros inteiros so criaturas do tipo -1, 0, e 32767. O tipo das strings que vamos usar chama-se unicode em Python 2.x, e str em Python 3.x. O tipo dos nmeros inteiros chama-se int em Python. A listagem 1.1 mostra o contedo do arquivo binint.py com alguns exemplos simples do uso de inteiros e strings. A legenda de cada listagem se inicia sempre com o nome do arquivo correspondente (aps um pouco de hesitao, eu decidi no disponibilizar estes arquivos; para testar os exemplos deste texto, voc precisar digitar os arquivos novamente: isto o (a) forar a praticar programao). A listagem um retrato el do arquivo, com duas excees: as palavras reservadas de Python esto sublinhadas na listagem (mas no no arquivo), e os espaos em branco dentro das strings esto enfatizados pelo smbolo . Ateno: os nmeros que aparecem esquerda da listagem no fazem parte do arquivo. Em Python, um comentrio inicia-se com #, e prossegue at o m de linha. A maior parte dos comandos de binint.py est explicada nos prprios comentrios. Alguns comentrios (sem inteno de trocadilho) adicionais, por linha: 1 Este comentrio especial torna o arquivo executvel em Linux. O caminho /usr/bin/python pode variar com a instalao, e mais ainda com o sistema operacional.

Matemtica Aplicada 2 3 4 8

14

Este comentrio informa que o arquivo est codicado em iso-8859-1, no bvio? Magia negra: tudo que escrito entre aspas passa a ser uma string to tipo unicode. Magia negra: print deixa de ser uma declarao, e passa a ser uma funo. Se voc no entendeu nada, no esquente. O operador ** signica exponenciao. Nesta linha, a varivel i torna-se um int, e recebe o valor 232 1. Este o maior valor sem sinal que um inteiro pode assumir, j que o maior inteiro que cabe em 32 bits 11111111111111111111111111111111 (binrio) = 1 231 + 1 230 + 1 229 + . . . + 1 22 + 1 21 + 1 20 = 4294967295 (decimal) .

1213 O comando for percorre em ordem um objeto itervel (no caso, a string aafro itervel, com a[0] == a , a[1] == , etc.). O corpo do for tem que ser indentado, e na listagem a indentao de 3 espaos. Desta forma, no caso da linha 13, o comando print( ord( ,c, ) = , ord(c)) executado 7 vezes, uma para cada caractere de a. A volta indentao anterior na linha 14 dene o m do for. Vamos agora sada do programa binint.py:aafro 7 4294967295 --- a : ord ( a ) = ord ( ) = ord ( a ) = ord ( f ) = ord ( r ) = ord ( ) = ord ( o ) = --- b : ord ( 4 ) = ord ( 2 ) = ord ( 9 ) = ord ( 4 ) = ord ( 9 ) = ord ( 6 ) = ord ( 7 ) = ord ( 2 ) = ord ( 9 ) = ord ( 5 ) =

97 231 97 102 114 227 111 52 50 57 52 57 54 55 50 57 53

print imprime com o formato apropriado inteiros e strings (e muito mais: quase tudo, de alguma forma!). O maior inteiro que cabe em 32 bits sem sinal 4294967295 (como j vimos acima). A posio do caractere a na tabela Unicode 97; a posio do caractere 231; a posio do caractere 4 52. Finalmente, o caractere Unicode de nmero 227 o . Se voc estiver vendo uma sada diferente da mostrada acima, com caracteres estranhos, no se assuste (demais): o seu arquivo binint.py e o terminal dentro do

15

1.2 Python Listagem 1.2: floats.py exemplo de uso de float e complex.

1 2 3 4 5 6 7 8 9 10 11 12 13 14

# !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n from math import pi , e , sin # pi , e , e seno # o seno de um nmero complexo vem com outro from cmath import sin as csin # nome , para no confundir from cmath import sqrt as csqrt # a raiz quadrada de um nmero complexo vem # com outro nome , para no confundir # imprime o valor de pi print ( pi = , pi ) print ( e = ,e ) # imprime o valor de e # imprime sen ( pi /2) print ( sen ( pi /2) = , sin ( pi /2)) i = csqrt ( -1.0) # neste programa , i == sqrt ( -1) print ( sen ( i ) = , csin ( i )) # imprime sen ( i )

qual voc est executando este programa certamente esto utilizando codicaes diferentes (veja se o apndice ?? pode ajudar). 1.2.2 Nmeros de ponto utuante reais e complexos Em primeiro lugar, um esclarecimento: no computador, no possvel representar todos os nmeros x R do conjunto dos reais, mas apenas um subconjunto dos racionais Q. Linguagens de programao mais antigas, como FORTRAN, ALGOL e PASCAL, mesmo assim chamavam estes tipos de REAL. A partir de C, e continuando com Python, em muitas linguagens passou-se a usar o nome mais adequado float. Python vem com uma grande quantidade de mdulos predenidos, e voc pode adicionar seus prprios mdulos. Deles, importam-se variveis e funes (e outras coisas) teis. Nosso primeiro exemplo do uso de nmeros de ponto utuante (float) e complexos (complex) no podia ser mais simples, na listagem 1.2 Eis a sada de floats.py:pi = 3.14159265359 e = 2.71828182846 sen ( pi /2) = 1.0 sen ( i ) = 1.17520119364 j

Os comentrios importantes seguem-se. Lembre-se de que na maioria dos casos voc est vendo aproximaes racionais de nmeros reais, e que o computador no pode lidar com todos os nmeros reais. Como era de se esperar, sen /2 = 1. Por outro lado, o seno de i um nmero puramente imaginrio, e vale 1,17520119364i. Note que Python usa a letra j para indicar um valor imaginrio (e no i). Na linha 13, eu forcei a barra, e criei a varivel que, em notao matemtica se escreveria i = 1. Mas eu tambm poderia ter simplesmente eliminado esta linha, e substitudo a linha 14 por print(sen(i) = ,csin(1j)) . Note que existem dois mdulos, math e cmath, para variveis reais e complexas, respectivamente. Em ambos, existe uma funo denominada sin, que calcula o seno. Para poder usar estas duas funes diferentes em meu programa floats.py, eu rebatizei a funo sin complexa de csin, no ato da importao do mdulo, com o mecanismo from ... import ... as ... (linha 6).

Matemtica Aplicada

16

Tabela 1.1: Vazes Mximas Anuais (m3 s1 ) no Rio dos Patos, PR, 19311999 Ano 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 1946 1947 1948 1949 1950 Vaz Mx 272.00 278.00 61.60 178.30 272.00 133.40 380.00 272.00 251.00 56.10 171.60 169.40 135.00 146.40 299.00 206.20 243.00 223.00 68.40 165.00 Ano 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 Vaz Mx 266.00 192.10 131.80 281.00 311.50 156.20 399.50 152.10 127.00 176.00 257.00 133.40 248.00 211.00 208.60 152.00 92.75 125.00 135.60 202.00 Ano 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 Vaz Mx 188.00 198.00 252.50 119.00 172.00 174.00 75.40 146.80 222.00 182.00 134.00 275.00 528.00 190.00 245.00 146.80 333.00 255.00 226.00 275.00 Ano 1991 1992 1993 1994 1995 1996 1997 1998 1999 Vaz Mx 131.00 660.00 333.00 128.00 472.00 196.00 247.50 451.00 486.00

Exerccios Propostos 1.3 Usando Python, converta 7777 da base 10 para a base 2. Sugesto: estude a documentao em www.python.org, e encontre a rotina pr-denida (built-in) que faz isto. 1.4 Como se faz para concatenar as strings bom e demais?

1.2.3 Obteno de uma curva de permanncia Uma funo distribuio acumulada (FDA) de probabilidade uma funo que nos informa qual a probabilidade de que uma varivel aleatria Q assuma um valor menor ou igual que um certo nivel q. Os valores de Q variam de experimento para experimento. Por exemplo, se Q a vazo mxima diria em um rio em um ano qualquer, o valor observado de Q varia de ano para ano. A tabela 1.1 d os valores da vazo mxima anual para o Rio dos Patos, PR, estao ANA (Agncia Nacional de guas do Brasil) 64620000, entre 1931 e 1999. Provavelmente, a maneira mais simples de se estimar uma FDA a partir de um conjunto de dados supor que os dados representam a totalidade das possibilidades, e que as observaes so equiprovveis (em analogia com os 6 nicos resultados possveis do lanamento de um dado no-viciado). No caso da tabela 1.1, se qi a vazo mxima do i-simo ano, teramos que a probabilidade de ocorrncia de qi P {Q = qi } = 1 , n (1.1)

17

1.2 Python

Listagem 1.3: patos-medmax.dat vazes mdia e mxima anuais, Rio dos Patos1931 1932 1933 1934 1935 21.57 25.65 4.76 11.46 28.10 272.00 278.00 61.60 178.30 272.00

onde n o nmero de observaes. Mas a FDA por denio F (q) = P {Q q}. (1.2) Para obt-la, preciso considerar os valores iguais ou menores que o valor de corte q. Portanto, ns devemos primeiro ordenar os qi s, de tal maneira que q0 q1 . . . qn1 . Note que a ordenao no altera (1.1). Aps a ordenao, o clculo de F (qi ) trivial: i 1 i+1 F (qi ) = = . (1.3) n n k=0 (1.3) chamada de distribuio acumulada emprica de probabilidade. Em Hidrologia, muitas vezes (1.3) denominada curva de permanncia. O resultado (i + 1)/n denominado uma posio de plotagem. Por diversos motivos, existem muitas outras posies de plotagem possveis para a FDA emprica. Uma muito popular (i + 1)/(n + 1). A discusso detalhada de posies de plotagem deve ser feita em um curso de Probabilidade e Estatstica, e no aqui, onde (1.3) serve (apenas) como um exemplo motivador. Os dados da tabela 1.1 esto digitados no arquivo patos-medmax.dat (Apndice A). Este arquivo contm 3 colunas contendo, respectivamente, o ano, a vazo mdia do ano, e a vazo mxima do ano. A listagem 1.3 mostra as 5 primeiras linhas do arquivo (que possui 69 linhas). O programa fqiemp.py, mostrado na listagem 1.4, calcula a curva de permanncia, ou FDA emprica, para as vazes mximas anuais do Rio dos Patos. Esta , simplesmente, uma tabela de duas colunas: a vazo observada (em ordem crescente), e o valor de (i + 1)/n. Como antes, os comentrios em fqiemp.py explicam muito do que est acontecendo. Mas h necessidade de explicaes adicionais: 5 Faz a diviso funcionar como em Python 3.x. Em Python 2.x, 3/4 == 0, e 3.0/4 == 0.75. Em Python 3.x, / sempre produz um resultado do tipo float (ou complex); e um novo operador // sempre produz diviso inteira. A funo pr-denida open abre o arquivo patos-medmax.dat, descrito acima e exemplicado na listagem 1.3. O segundo argumento de open, a string rt , diz duas coisas: com r, que se trata de uma operao de leitura, (read), de um arquivo que j existe: isto garante que patos-medmax.dat no ser modicado por fqiemp.py; com t, que se trata de um arquivo texto. Um arquivo texto um arquivo formado por linhas, sendo cada linha uma string. As linhas so separadas por caracteres (invisveis no editor de texto) de m de linha, que convencionalmente ns indicamos por \n. Um arquivo texto um objeto itervel, que pode ser acessado linha a linha.

6

Matemtica Aplicada

18

Listagem 1.4: fqiemp.py clculo de uma FDA emprica1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n from __future__ import division fin = open ( patos - medmax . dat , rt ) # qmax = [] # # for linha in fin : campo = linha . split () # # print ( campo ) qi = float ( campo [2]) # qmax . append ( qi ) # fin . close () # fou = open ( fqiemp . dat , wt ) # qmax . sort () # n = len ( qmax ) # # for i in range ( n ): qi = qmax [ i ] # Fi = ( i +1)/ n # fou . write ( %8.2 f %8.6 f \ n % ( qi , Fi )) # fou . close () #

abre o arquivo de dados ( entrada ) uma lista vazia loop nas linhas do arquivo separa os campos para ver campo a campo na tela a vazo o terceiro campo adiciona qi lista qmax fecha o arquivo de entrada abre o arquivo de sada ordena a lista tamanho da lista loop nos elementos da lista ordenada vazo posio de plotagem imprime uma linha fim de papo

7

Nesta linha, qmax inicializado como uma lista vazia. Uma lista uma sequncia de objetos quaisquer. Dada uma lista a, seus elementos so a[0], a[1], etc.. Uma lista a com n objetos vai de a[0] at a[n-1]. Este o loop de leitura do programa. linha uma string que contm em cada iterao uma das linhas de patos-medmax.dat. preciso separar uma linha (veja a listagem 1.3) em seus 3 campos. O mtodo1 split separa a string linha em 3 strings, e as coloca em uma lista campo == [ campo[0], campo[1], campo[2] ], usando os espaos em branco como separadores (que so eliminados). Cada um dos campos agora ele mesmo uma string, com os valores separados. As 5 primeiras linhas que aparecem na tela devido ao comando print da linha 10 so:[ [ [ [ [ 1931 1932 1933 1934 1935 , , , , , 21.57 , 272.00 25.65 , 278.00 4.76 , 61.60 ] 11.46 , 178.30 28.10 , 272.00 ] ] ] ]

8 9

10

11 12 13 141

No entato, estes campos ainda so strings, e no floats. Aqui o terceiro campo convertido em um float e vai para a varivel qi. Finalmente, qi includa na lista qmax. bom estilo de programao fechar cada arquivo previamente aberto. Agora o programa abre o arquivo de sada, fqiemp.dat. w signica abertura para escrita (write); e t que o arquivo ser texto.

Em Python, um mtodo uma rotina que pertence a uma varivel ou um tipo (a rigor, a uma classe, mas este no um curso de programao)

19

1.2 Python

Listagem 1.5: fqiemp.dat FDA emprica da vazo mxima anual no Rio dos Patos.56.10 61.60 68.40 75.40 92.75 0.014493 0.028986 0.043478 0.057971 0.072464

15 16 17 18 19 20

sort um mtodo pr-denido para qualquer lista, que a ordena (por default em ordem crescente). len uma funo pr-denida que retorna o nmero de elementos da lista. loop para impresso no arquivo de sada. Obtm o i-simo elemento de qmax, colocado na varivel qi, que re-utilizada para este m. Calcula a posio de plotagem, utilizando o operador / para dividir dois ints e gerar um resultado correto do tipo float (veja comentrio sobre a linha 5). write um mtodo do arquivo fou. write sempre escreve seu nico argumento, que tem que ser uma string, no arquivo. Trata-se portanto do problema inverso do loop de leitura, que transformava strings em floats: agora precisamos transformar floats em uma string. para isto que serve o operador %: ele tem sua esquerda uma string com os campos de formatao especiais %8.2f e %8.6f; sua direita uma tupla 2 (qi,Fi) com tantos elementos quantos so os campos de formatao. O primeiro elemento ser substitudo no primeiro campo de formatao por uma string com 8 caracteres, sendo 2 deles para casas decimais. O segundo elemento ser substitudo no segundo campo de formatao por uma string com 8 caracteres, sendo 6 deles para casas decimais. A string resultante, que precisa conter explicitamente o caractere de m de linha \n, ser escrita no arquivo de sada.

As primeiras 5 linhas do arquivo de sada fqiemp.dat so mostradas na listagem 1.5.Exerccios Propostos 1.5 Dado um nmero inteiro p lido do terminal, escreva um programa Python para procurar o ndice i de uma lista a de 10 nmeros inteiros denida internamente no programa tal que a[i] == p, e imprimir o resultado.

1.2.4 Arquivos texto e arquivos binrios Arquivos texto so legveis por seres humanos. Qualquer representao interna primeiramente traduzida para uma string de caracteres antes de ser escrita em um arquivo texto. Arquivos binrios em geral armazenam informao com a mesma representao interna utilizada pelo computador para fazer contas, etc..2

Em Python uma tupla uma lista imutvel

Matemtica Aplicada

20

Listagem 1.6: bintext.py Exemplo de arquivo texto e arquivo binrio1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import print _functio n i = 673451 bdet = bin ( i ) print ( bdet ) fot = open ( i . txt , wt ) fot . write ( %6 d % i ) fot . close () from struct import pack b = pack ( i ,i ) fob = open ( i . bin , wb ) fob . write ( b ) fob . close () print ( type ( i )) print ( type ( b )) from os . path import getsize print ( tamanho do arquivo txt = % d bytes print ( tamanho do arquivo bin = % d bytes

% getsize ( i . txt )) % getsize ( i . bin ))

Como sempre, um exemplo vale por mil palavras. Na listagem 1.6, temos o programa bintext.py, que produz dois arquivos: o arquivo texto i.txt, e o arquivo binrio i.bin. Cada um destes arquivos contm o nmero inteiro i == 673451. Eis aqui a dissecao de bintext.py: 5 8 10 A funo pr-denida bin produz a representao do objeto na base 2. Escreve a varivel i no arquivo i.txt, usando um campo de 6 caracteres, que o tamanho necessrio para escrever i na base 10. A nica maneira de ler ou escrever um arquivo binrio em Python por meio de byte strings: strings em que cada elemento um byte (8 bits). O mdulo struct prov a converso de, e para, byte strings. Converte i em uma byte string. Como a representao interna de i em 4 bytes, o tamanho de b ser de 4 bytes.

11

1516 Verica os tipos de i e b. 1819 Mostra o tamanho de cada um dos arquivos gerados. Finalmente, a sada de bintext.py 0 b10100100011010101011 < type int > < type str > tamanho do arquivo txt = 6 bytes tamanho do arquivo bin = 4 bytes

Conte os 32 bits na primeira linha (aps o prexo 0b, que indica a representao na base 2): eles correspondem aos 4 bytes que custa ao programa para guardar o nmreo inteiro 673451. Repare que o arquivo binrio menor que o arquivo texto. Em geral, arquivos binrios tendem a ser menores (para a mesma quantidade de informao). A outra grande vantagem que a leitura e a escritura de arquivos binrios muito mais rpida, porque no h necessidade de traduzir a informao de, e para, strings.

21

1.2 Python

Listagem 1.7: writearr.py Escreve um arquivo binrio contendo 3 linhas, cada uma das quais com um array de 10 floats.1 2 3 4 5 6 7 8 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from numpy . random import rand fob = open ( a . bin , wb ) for k in range (3): a = rand (10) a . tofile ( fob ) fob . close ()

# # # #

abre um loop em gera um escreve

arq binrio para escrita 3 " linhas " array com 10 nmeros aleatrios uma " linha "

Na prtica, arquivos binrios esto invariavelmente associados ao uso de arrays, pelo menos em Engenharia. O mdulo Numpy3 , que no vem com Python, e necessita ser instalado parte, proporciona um tipo chamado array, e permite manipular com boa ecincia computacional vetores e matrizes em Python. O tipo permite na prtica substituir listas (tipo list); ademais, tudo ou quase tudo que funciona com listas tambm funciona com arrays. Neste texto, ns faremos a partir de agora amplo uso de Numpy e de seu tipo array. A referncia completa de Numpy est disponvel em domnio pblico (Oliphant, 2006), podendo tambm ser comprada pela Internet. Alm de denir arrays, Numpy tambm proporciona seus prprios mtodos e funes para ler e escrever de e para arquivos binrios. Vamos ento dar dois exemplos muito simples, porm muito esclarecedores. Primeiro, um programa para escrever um array. A listagem 1.7 mostra o programa writearr.py. O programa usa uma rotina disponvel em numpy, rand, para devolver um array com 10 nmeros aleatrios entre 0 e 1. writearr.py repete por 3 vezes a gerao de a e a sua escritura no arquivo binrio a.bin: a escritura utiliza o mtodo tofile. Repare que tofile um mtodo do array a; ele no precisa ser importado, pois ele j faz parte da varivel a a partir do momento em que ela criada. writearr roda silenciosamente: no h nenhuma sada na tela. No entanto, se procurarmos no disco o arquivo gerado, teremos algo do tipo> ls -l a . bin -rw -r - -r - - 1 nldias nldias 240 2011 -08 -28 14:08 a . bin

O arquivo gerado, a.bin, possui 240 bytes. Em cada uma das 3 iteraes de writearr.py, ele escreve 10 floats no arquivo. Cada float custa 8 bytes, de modo que em cada iterao 80 bytes so escritos. No nal, so 240. importante observar que o arquivo a.bin no possui estrutura: ele no sabe que dentro dele mora o array a; ele , apenas, uma linguia de 240 bytes. Cabe a voc, programadora ou programador, interpretar, ler e escrever corretamente o arquivo. Prosseguimos agora para ler o arquivo binrio gerado. Isto feito com o programa readarray.py, mostrado na listagem 1.8. O programa importa a rotina fromfile de Numpy, a partir da qual 3 instncias de a so lidas do arquivo a.bin e impressas com formato na tela. Eis a sua sada:0.2327 0.6117 0.7713 0.5942 0.1799 0.3156 0.1473 0.4299 0.0870 0.0846 0.4301 0.9779 0.0322 0.4833 0.6097 0.4387 0.0639 0.1399 0.4350 0.7737 0.5809 0.0382 0.6567 0.8062 0.8427 0.2511 0.2897 0.5785 0.2892 0.03853

numpy.scipy.org

Matemtica Aplicada

22

Listagem 1.8: readarr.py Escreve um arquivo binrio contendo 3 linhas, cada uma das quais com um array de 10 floats.1 2 3 4 5 6 7 8 9 10 11 12 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import print _functio n from numpy import fromfile # para escrever na tela from sys import stdout fib = open ( a . bin , rb ) # abre o arquivo binrio for k in range (3): # loop em 3 " linhas " a = fromfile ( fib , float ,10) # l um array com 10 floats # imprime com formato na tela for e in a : stdout . write ( %5.4 f % e ) stdout . write ( \ n ) fib . close ()

O que vemos so os nmeros aleatrios das 3 instncias de a escritas pelo programa writearr.py.

1.3 GnuplotMas como a cara da funo distribuio acumulada emprica? Ns gostaramos de ver a funo, isto , plot-la. para isto que serve Gnuplot. Um programa Gnuplot especica o arquivo de onde vm os dados; que colunas queremos utilizar; que smbolos ou tipos de linhas, bem como espessuras ou tamanhos ns desejamos, e o nome do arquivo de sada (na falta de um nome para o arquivo de sada, uma tela com o grco dever se abrir para visualizao interativa). Via de regra, o arquivo de sada de um programa Gnuplot no um arquivo texto, mas um arquivo com uma descrio grca qualquer. Extenses comuns de nomes de arquivos que indicam imagens so .jpg, .png, .eps, e .emf (esta ltima em Windows). Gnuplot permite gerar arquivos de imagem em qualquer um destes formatos (e muitos outros!). Neste texto, ns vamos sempre gerar pares de arquivos (quase idnticos) .eps e .emf. Em Linux, um arquivo a.eps facilmente conversvel no arquivo a.pdf mediante o comando epstopdf a.eps ; em Windows, um arquivo a.emf est pronto para ser visualizado (com um clique), inserido em documentos, etc.. Ns j temos a FDA emprica da vazo mxima anual no Rio dos Patos, no arquivo fqiemp.dat; agora, ns vamos plot-la com o programa Gnuplot fiqemp.plt: Como sempre, preciso explicar: 5 6 7 8 9 Gera um arquivo de sada no padro Encapsulated Postscript, em preto-ebranco, com tipo Times Roman, no tamanho 18pt. O nome do arquivo de sada. Gera um arquivo que mostra as convenes utilizadas neste terminal. Este arquivo pode ser visto na gura 1.1 O nome do arquivo de sada. Formato do eixo y, que usa a mesma conveno de Python (a bem da verdade, de C).

23

1.3 Gnuplot

Listagem 1.9: fiqemp.plt programa para plotar a FDA emprica do Rio dos Patos1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 set encoding iso_8859_1 # ------------------------------------------------------------------------------# no mundo de linux # ------------------------------------------------------------------------------set terminal postscript eps monochrome Times - Roman 18 set output test . eps test set output fqiemp . eps set format y %3.1 f set grid set xlabel Vazo mxima anual ( m3 / s ) set ylabel FDA emprica plot fqiemp . dat using 1:2 notitle with points pt 7 # -----------------------------------------------------------------------------# agora no mundo de Windows # -----------------------------------------------------------------------------set terminal emf monochrome Times New Roman 18 set output test . emf test set output fqiemp . emf replot exit

10 11 12 13 17 18 19 20 21

Mostra o reticulado da escala. Nome do eixo x. Nome do eixo y. Nome do arquivo de dados a ser utilizado; colunas para x e y; use pontos, com smbolo 7 (veja a conveno na gura 1.1): gera o arquivo fqiemp.eps. Mude o tipo de arquivo de sada para emf, usando o tipo Times New Roman, no tamanho 18 pt. O nome do arquivo de sada. O mesmo que a linha 7. O nome do arquivo de sada. Plota novamente, usando todas as denies que no tiverem sido mudadas (s mudamos o terminal de sada); agora, gera o arquivo fqiemp.emf.

Finalmente, podemos ver a FDA emprica gerada por fqiemp.plt, na gura 1.2. Uma ltima observao de cunho esttico. Neste texto, as guras que so geradas por programas Gnuplot descritos explicitamente, e cujas listagens .plt so dadas no texto, utilizam o tipo postscript Times-Roman. Olhe bem para o tipo utilizado nas guras 1.1 e 1.2: este no o tipo utilizado no restante do texto (por exemplo, compare com as legendas). J as guras geradas para m ilustrativo utilizam o mesmo tipo que o texto; veja por exemplo as guras 2.1, e ??. Utilizar um tipo de letra diferente nas guras e no texto ligeiramente inconsistente, do ponto de vista tipogrco. Entretanto, a escolha de Times-Roman para

Matemtica AplicadaTerminal Test show ticscale 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 9 22 23 24

24

(color) filled polygon:

rotated ce+ntred text

te d

by

+4

left justified centre+d text right justified test of character width: 12345678901234567890

ro

ta

linewidth lw 6 lw 5 lw 4 lw 3 lw 2 lw 1 pattern fill 3 4 5 6

5 ro ta te d by 4 5 de g

de g

0

1

2

7

8

Figura 1.1: Convenes do terminal postscript eps monochrome 18

Times-Roman

a sada em postscript, e de Times New Roman para a sada em EMF, foi feita em nome da simplicidade e da universalidade. praticamente impossvel no existir o tipo Times-Roman em seu sistema operacional Linux, e praticamente impossvel no existir o tipo Times New Roman em seu sistema operacional Windows. Na preparao de um texto mais bem-elaborado tipogracamente voc deve, se possvel, utilizar o mesmo tipo nas guras e no texto. Mas isto, alm de ser outra histria, do gosto do fregus. 1.3.1 Ajuste de uma funo A denio da FDA emprica em (1.3) muito limitada: bvio que nenhum destes 69 valores da tabela 1.1, ao contrrio do que acontece no jogo de dados, jamais se repetir. Tambm bvio que valores maiores do que mximo de 660 m3 s1 podero ocorrer no futuro (em relao a 1999 . . . ), e isto uma poderosa fonte de crtica posio de plotagem ingnua que ns utilizamos em (1.3): ela prev que a probabilidade de uma vazo mxima anual maior do que 660 m3 s1 nula! Em suma, nossa esperana que uma relao analtica, uma frmula, faa um trabalho no mnimo um pouco melhor do que a FDA emprica foi capaz de fazer para prever as probabilidades da vazo mxima anual no Rio dos Patos mas espere para ver o resultado! Vamos tentar ajustar uma frmula aos dados. A frmula que escolheremos a da distribuio de Weibull, F (q) = 1 exp (q/)k . (1.4)

No existe absolutamente nenhum motivo transcedental para utilizar (1.4): ela apenas atende aos requisitos fundamentais de uma FDA (tende a 0 em x ; e a 1 em x +). A rigor, a Weibull denida apenas para q 0, que o que se espera de qualquer vazo em rios nenhum rio tem vazo negativa; nenhum rio corre do mar para suas cabeceiras.

25

1.3 Gnuplot

1.0 0.9 0.8 0.7 FDA emprica 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0 100 200 300 400 500 Vazo mxima anual (m3/s) 600 700

Figura 1.2: FDA da vazo mxima anual, Rio dos Patos

Listagem 1.10: weibull.plt como ajustar uma FDA analtica (Weibull) aos dados da FDA emprica1 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 set encoding iso_8859_1 # ------------------------------------------------------------------------------# no mundo de linux # ------------------------------------------------------------------------------set terminal postscript eps monochrome Times - Roman 18 # ------------------------------------------------------------------------------# ajusta lambda e k # ------------------------------------------------------------------------------lambda = 1000.0 k = 1.0 F ( x ) = 1.0 - exp ( -(( x / lambda )** k )) fit F ( x ) fqiemp . dat using 1:2 via lambda , k set output weibullfit . eps set format y %3.1 f set grid set xlabel Vazo mxima anual ( m3 / s ) set ylabel FDA set xrange [0:1000] set yrange [0:1] set xtics 0 ,100 set ytics 0 ,0.1 set samples 10000 plot fqiemp . dat using 1:2 notitle with points pt 7 , \ F ( x ) with lines notitle lt 1 lw 3 # -----------------------------------------------------------------------------# agora no mundo de Windows # -----------------------------------------------------------------------------set terminal emf monochrome Times New Roman 18 set output weibullfit . emf replot exit

Matemtica Aplicada

26

No programa weibullfit.plt, mostrado na listagem 1.10, veja como Gnuplot lida facilmente com o ajuste dos valores de e de k aos dados de que dispomos. Gnuplot produz informaes na tela sobre o ajuste requisitado pelo comando fit; as ltimas linhas da sada em tela soAfter 11 iterations the fit converged . final sum of squares of residuals : 0.0604353 rel . change during last iteration : -6.18025 e -06 degrees of freedom ( FIT_NDF ) rms of residuals ( FIT_STDFIT ) = sqrt ( WSSR / ndf ) variance of residuals ( reduced chisquare ) = WSSR / ndf Final set of parameters ======================= lambda k = 233.486 = 2.62023 : 67 : 0.0300336 : 0.000902019

Asymptotic Standard Error ========================== +/ - 1.283 +/ - 0.05594 (0.5495%) (2.135%)

correlation matrix of the fit parameters : lambda k lambda k 1.000 -0.346 1.000

donde = 233,486, k = 2,62023. Na listagem 1.10, comentrios adicionais so: 910 Os valores iniciais de lambda e k inuenciam nos valores nais encontrados ou no. Por exemplo, se mudarmos a linha 9 para lambda = 1.0 , o comando fit da linha 12 no consegue ajustar a curva F (x); muitas vezes voc precisar fazer vrias tentativas com os valores iniciais at ter sucesso. 11 22 23 possvel denir funes que sero depois plotadas. possvel aumentar o nmero de pontos utilizados para traar curvas, quando a aparncia das mesmas for quebrada. possvel aumentar a espessura das curvas, com lw (line width): o nmero a seguir especica a espessura desejada: veja a gura 1.1.

A gura 1.3 mostra o grco gerado: note como a cauda da direita, ou seja, os valores maiores de vazo mxima, so mal ajustados: F (q) subestima o valor de q para probabilidades altas (prximas de 1): o uso da Weibull para este conjunto de dados seria desastroso, uma vez que as maiores cheias seriam subestimadas!Exerccios Propostos 1.6 Usando Gnuplot, plote em tons de cinza a rea debaixo da curva y = x para x = 2. Quanto vale esta rea?

1.4 Python de novo: projeto probabilsticoDe fato, um projeto probabilstico a denio de um valor de corte q cuja probabilidade de excedncia seja menor que um certo risco predeterminado . Suponha que, para o Rio dos Patos, ns desejemos calcular o valor q da vazo mxima anual

271.0 0.9 0.8 0.7 0.6 FDA 0.5 0.4 0.3 0.2 0.1 0.0 0 100 200 300

1.4 Python de novo: projeto probabilstico

400 500 600 700 Vazo mxima anual (m3/s)

800

900

1000

Figura 1.3: Ajuste de uma FDA Weibull aos dados de vazo mxima anual do Rio dos Patos cuja probabilidade de excedncia em um ano qualquer seja = 0.05 (5%). Isto signica que precisamos resolver a equao P {Q > q} = 1 P {Q q} = , 1 F (q) = , F (q) = 1 ,

(1.5)

para q. Quando a F (q) a FDA emprica, no existe uma equao para se resolver, mas sim uma tabela na qual devemos interpolar o valor desejado de q. Primeiramente ns precisamos escrever uma funo que, dado um valor xc, procure o intervalo em que ele se encontra na coluna x da tabela, e em seguida calcule o yc correspondente interpolando linearmente entre os valores da coluna y para o mesmo intervalo. O mdulo interp, arquivo interp.py, mostrado na listagem 1.11, faz isto. interp no um programa: ele um mdulo que contm uma nica funo (que possui tambm o nome interp). Esta rotina estar disponvel para qualquer programa Python por meio de uma declarao do tipo from interp import interp . A dissecao do mdulo e da funo interp 6 9 Em Python, uma funo denida com a palavra reservada def. assert signica assegure-se de que; == um operador lgico que retorna True se os operandos forem iguais.

1718 iu para upper, il para lower: so os ndices do intervalo que contm xc. O intervalo se inicia igual ao intervalo da tabela toda.

Matemtica Aplicada

28

Listagem 1.11: interp.py mdulo com funo para calcular uma interpolao linear.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 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n from __future__ import division def interp ( xc , x , y ): nx = len ( x ) ny = len ( y ) assert nx == ny # x e y devem ser listas com tamanhos iguais # deve haver pelo menos 2 pontos if nx < 2: exit ( " interp : n = % d deve ser >= 2\ n " % nx ) if ( xc < x [0]) or ( xc > x [ nx -1]): # xc deve estar no intervalo da lista x exit ( " xc = % lf fora de alcance : % f % f \ n " % ( xc , x [0] , x [ nx -1])) # -----------------------------------------------------------------------------# inicia uma busca binria # -----------------------------------------------------------------------------iu = nx - 1 il = 0 while ( iu - il > 1): im = ( iu + il )//2 if ( xc 0) $ assume ( lam > 0) $ declare ( [ k ] , noninteger ) $ F : 1 - exp ( -( x / lam )^ k ) $ fdp : diff (F , x ) $ mu : integrate ( x * fdp ,x ,0 , inf ) ;

6

Dene F (F (x)): note que no F(x)! Note tambm que a atribuio, em Maxima, no feita com = (como em Python ou Gnuplot), mas sim com :; = ca reservado para igualdade. armazena a derivada de F (x), que em probabilidade se chama funo densidade de probabilidade, na varivel fdp. Calcula a integral denidora da mdia. O resultado de rodar mulambdak.max, com o comando maxima -b mulambdak.max

7 8

Maxima 5.21.1 http :// maxima . sourceforge . net using Lisp GNU Common Lisp ( GCL ) GCL 2.6.7 ( a . k . a . GCL ) Distributed under the GNU Public License . See the file COPYING . Dedicated to the memory of William Schelter . The function bug_report () provides bug reporting information . (% i1 ) batch ( mulambdak . max ) read and interpret file : # p / home / nldias / work / graduacao / matap / aulas / mulambdak . max (% i2 ) declare ([ k ] , real ) (% i3 ) declare ([ lam ] , real ) (% i4 ) assume ( k > 0) (% i5 ) assume ( lam > 0) (% i6 ) declare ([ k ] , noninteger ) x k (% i7 ) F : 1 - exp ( - ( - - -) ) lam (% i8 ) fdp : diff (F , x ) (% i9 ) mu : integrate ( x fdp , x , 0 , inf ) k + 1 (% o9 ) gamma ( - - - - -) lam k (% o9 ) mulambdak . max

Na linha (%o9), gamma signica a funo gama (x), que ns vamos encontrar muitas vezes neste curso, mas que no vem ao caso detalhar agora. O resultado analtico que ns obtivemos com Maxima = Exerccios Propostos 1.9 Com Maxima, calcule a derivada de f (x) = ln(sen(ex )).

k+1 . k

(1.8)

Matemtica Aplicada1.10 Com Maxima, calcule1

32

1 x3/2

dx.

1.11 Com Maxima, obtenha as razes de 3 29 x3 + x2 x + 15 = 0. 2 2

2Um pouco de polinmios, integrais, sries . . .2.1 Integrao numrica: motivaoSuponha que voc deseje traar uma curva com as seguintes propriedades: passar pelos pontos (0, 0), (1, 5/2), (2, 7/2) e (4, 4); possuir derivada igual a 0 em x = 4. Existem muitas curvas com estas propriedades, mas uma candidata natural um polinmio de grau 5, uma vez que existem 5 propriedades ou graus de liberdade na lista acima. Portanto, f (x) = ax4 + bx3 + cx2 + dx + e, df , g(x) = dx f (0) = 0, f (1) = 5/2, f (2) = 7/2, f (4) = 4, g(4) = 0. Agora, podemos obter facilmente a, b, c, d, e com Maxima, com o programa achapol.max da listagem 2.1. A dissecao de achapol.max a seguinte: 12 37 8 Dene f (x) e g(x). Dene o valor de f (x) em x = 0, 1, 2, 4, e o valor de g(x) em x = 4, em funo de a, b, c, d, e. Resolve um sistema linear 5 5 em a, b, c, d, e. A sada de achapol.max mostrada na listagem 2.2 Portanto, 1 13 17 11 f (x) = x4 + x3 x2 + x. 48 48 12 3 33

(2.1)

Matemtica Aplicada

34

Listagem 2.1: achapol.max Polinmio com propriedades denidas1 2 3 4 5 6 7 8 f : a * x ^4 + b * x ^3 + c * x ^2 + d * x + e ; g : diff (f , x ); eq1 : f , x =0 ; eq2 : f , x =1 ; eq3 : f , x =2 ; eq4 : f , x =4 ; eq5 : g , x =4 ; solve ([ eq1 = 0 , eq2 = 5/2 , eq3 = 7/2 , eq4 = 4 , eq5 = 0] , [a ,b ,c ,d , e ]) ;

Listagem 2.2: Sada de achapol.maxMaxima 5.21.1 http :// maxima . sourceforge . net using Lisp GNU Common Lisp ( GCL ) GCL 2.6.7 ( a . k . a . GCL ) Distributed under the GNU Public License . See the file COPYING . Dedicated to the memory of William Schelter . The function bug_report () provides bug reporting information . (% i1 ) batch ( achapol . max ) read and interpret file : # p / home / nldias / work / graduacao / matap / aulas / achapol . max 2 3 4 (% i2 ) f : e + d x + c x + b x + a x 4 3 2 (% o2 ) a x + b x + c x + d x + e (% i3 ) g : diff (f , x ) 3 2 (% o3 ) 4 a x + 3 b x + 2 c x + d (% i4 ) ev ( eq1 : f , x = 0) (% o4 ) e (% i5 ) ev ( eq2 : f , x = 1) (% o5 ) e + d + c + b + a (% i6 ) ev ( eq3 : f , x = 2) (% o6 ) e + 2 d + 4 c + 8 b + 16 a (% i7 ) ev ( eq4 : f , x = 4) (% o7 ) e + 4 d + 16 c + 64 b + 256 a (% i8 ) ev ( eq5 : g , x = 4) (% o8 ) d + 8 c + 48 b + 256 a 5 7 (% i9 ) solve ([ eq1 = 0 , eq2 = -, eq3 = -, eq4 = 4 , eq5 = 0] , [a , b , c , d , e ]) 2 2 1 13 17 11 (% o9 ) [[ a = - --, b = --, c = - --, d = --, e = 0]] 48 48 12 3 (% o9 ) achapol . max

Suponha agora que voc deseje calcular I=5 1

f (x) dx.

claro que esta integral pode ser calculada analiticamente com Maxima:(% i1 ) [ a : -1/48 , b : 13/48 , c : -17/12 , d : 11/3] ; 1 13 17 11 (% o1 ) [ - --, --, - --, - -] 48 48 12 3 (% i2 ) f : a * x ^4 + b * x ^3 + c * x ^2 + d * x ; 4 3 2 x 13 x 17 x 11 x (% o2 ) - -- + ----- - ----- + ---48 48 12 3 (% i3 ) integrate (f ,x ,1 ,5) ; 1321 (% o3 ) ---90

35f (x)

2.1 Integrao numrica: motivao

4 f (x)

3

2

h(x)

1

0

0

1

2

3

4

5

x

Figura 2.1: Integrao numrica de uma funo(% i4 ) bfloat (%); (% o4 )

1 . 4 6 7 7 7 7 7 7 7 7 7 7 7 7 8 b1

Logo, I = 1321/90 14,6778 (ateno para o arredondamento). Entretanto, nem todas as integrais podem ser calculadas assim, por uma srie de motivos: a funo pode no possuir uma integral em forma fechada, a funo pode ser denida por pontos, mas no por uma frmula, a funo pode ser difcil de integrar analiticamente, e voc pode querer ter uma idia inicial do valor da integral, etc.. Suponha portanto que voc no soubesse que I = 1321/90, mas fosse capaz de calcular em princpio tantos pontos de f (x) quantos fossem necessrios: qual a sua aposta para o valor de I? Considere a gura 2.1: uma das aproximaes mais bvias mas tambm mais grosseiras substituir a funo por uma reta ligando os pontos (1, f (1)) e (5, f (5)). A rea do trapzio a nossa primeira aproximao para a integral: I0 = f (1) + f (5) 4 = 12,50. 2

Nada mal, considerando o valor verdadeiro 14,6778! Mas pode ser melhorada, com o uso de dois trapzios. Para isto, basta calcular f (3) e somar as reas dos dois trapzios resultantes: I1 = f (1) + f (3) f (3) + f (5) 2+ 2 = 14,000. 2 2

Note que estamos muito prximos do valor verdadeiro com apenas 2 trapzios. Mas existe uma alternativa analtica mais inteligente do que trapzios: como temos 3 pontos, ns podemos aproximar a curva por uma parbola do 2 o grau passando

Matemtica Aplicada

36

por estes 3 pontos. O programa passaquad.max, mostrado na listagem 2.3, faz este trabalho, e acha os coecientes a, b, c da parbola h(x) = ax2 + bx + c que passa por (1, f (1)), (3, f (3)) e (5, f (5)). No embalo, passaquad.max redene h(x) com os coecientes encontrados, e j calcula a integral de h(x) entre 1 e 5. A parbola h tambm mostrada na gura 2.1. Listagem 2.3: passaquad.max parbola h(x) = ax2 +bx+c passando por (1, f (1)), (3, f (3)) e (5, f (5)).1 2 3 4 5 6 7 8 9 10 11 12 f : ( -1/48)* x ^4 + (13/48)* x ^3 - (17/12)* x ^2 + (11/3)* x ; y1 : f , x =1 $ y3 : f , x =3 $ y5 : f , x =5 $ h : a * x ^2 + b * x + c$ eq1 : ev (h , x =1) = y1 ; eq2 : ev (h , x =3) = y3 ; eq3 : ev (h , x =5) = y5 ; solve ( [ eq1 , eq2 , eq3 ] ,[a ,b , c ]) ; abc : % ; h : h , abc ; integrate (h ,x ,1 ,5);

A sada de passaquad.max mostrada na listagem 2.4. Com os coecientes a, b, c de h(x), ns podemos calcular analiticamente nossa prxima aproximao: I2 = =2 1 1

h(x) dx2

29 = = 14,5000. 2

3 2 23 5 x + x+ dx 16 16 4

At agora ns avaliamos 3 alternativas de integrao numrica de f (x): com um trapzio, com dois trapzios, e com uma parbola. A tabela 2.1 d um resumo dos resultados alcanados. O erro relativo de cada estimativa = Ik I . I (2.2)

Uma nica parbola foi capaz de estimar I com um erro relativo ligeiramente superior a 1%. Um caminho geral para a integrao numrica est aberto: aumentar o nmero de elementos de integrao (no nosso caso foram trapzios) e/ou aumentar a ordem do polinmio aproximador da funo por um certo nmero de pontos. Este o contedo da prxima seo.Exerccios Propostos 2.1 Se f (x) = sen x, qual a estimativa de I = 0

f (x) dx = 2 com um trapzio ?

2.2 Provavelmente, voc no est muito contente com o resultado do Problema 2.1.

37

2.2 A regra do trapzio Listagem 2.4: Sada de passaquad.max

batching / home / nldias / work / graduacao / matap / aulas / passaquad . max 2 3 4 11 x 17 x 13 x ( - 1) x (% i2 ) f : ---- - ----- + ----- + - - - - - - - 3 12 48 48 4 3 2 x 13 x 17 x 11 x (% o2 ) - -- + ----- - ----- + ---48 48 12 3 (% i3 ) ev ( y1 : f , x = 1) (% i4 ) ev ( y3 : f , x = 3) (% i5 ) ev ( y5 : f , x = 5) 2 (% i6 ) h : c + b x + a x (% i7 ) eq1 : ev (h , x = 1) = y1 5 (% o7 ) c + b + a = 2 (% i8 ) eq2 : ev (h , x = 3) = y3 31 (% o8 ) c + 3 b + 9 a = -8 (% i9 ) eq3 : ev (h , x = 5) = y5 15 (% o9 ) c + 5 b + 25 a = -4 (% i10 ) solve ([ eq1 , eq2 , eq3 ] , [a , b , c ]) 3 23 5 (% o10 ) [[ a = - --, b = --, c = -]] 16 16 4 (% i11 ) abc : % 3 23 5 (% o11 ) [[ a = - --, b = --, c = -]] 16 16 4 (% i12 ) ev ( h : h , abc ) 2 3 x 23 x 5 (% o12 ) - ---- + ---- + 16 16 4 (% i13 ) integrate (h , x , 1 , 5) 29 (% o13 ) -2 (% o13 ) passaquad . max

a) Aproxime a integral I do Problema 2.1 com dois trapzios, entre x = 0 e /2, e entre x = pi/2 e . b) Aproxime a integral pela integral da parbola quadrtica g(x) = ax2 + bx + c passando pelos pontos (0, 0), (/2, 1) e (, 0). c) Aproxime a integral de f (x) pela integral da parbola cbica h(x) = ax3 + bx2 + cx + d passando pelos mesmos pontos acima, e com h (/2) = 0. 2.3 Com Gnuplot, plote f (x) = sen x e a aproximao g(x) obtida no Problema 2.2.

2.2 A regra do trapzioVamos comear a melhorar nossas estimativas de I pelo mtodo de fora bruta, de aumentar o nmero de trapzios. Isto nos levar ao mtodo talvez mais simples de integrao numrica que vale a pena mencionar, denominado Regra do Trapzio.

Matemtica Aplicada Tabela 2.1: Estimativas numricas de I e seus erros relativos . Integral Exato Um trapzio Dois trapzios Uma parbolaf (x)

38

Valor

14,6778 0 12,5000 0,1483 14,0000 0,0461 14,5000 0,0121

f (x)

x0

x1

x2

x3

x4

x5

x6

x7

x8

x

Figura 2.2: A regra do trapzio, com n = 4 e n = 8 trapzios. A gura 2.2 mostra a mesma funo f (x) da seo anterior; agora, entretanto, ns desenhamos 4 e 8 trapzios sob a curva f (x). evidente que a rea sob f (x) est agora muito bem aproximada com n = 8 trapzios. O seu valor pode ser estimado por I3 = com x0 = 1, xn = 5, xn x0 x = . n Prosseguindo, (2.3) pode ser re-escrita: I3 = f (x0 ) + f (x1 ) f (x1 ) + f (x2 ) f (xn1 ) + f (xn ) + + ... + x 2 2 2 x = f (x0 ) + 2 f (x1 ) + f (x2 ) + . . . + f (xn1 ) + f (xn ) 2 x = (Se + 2Si ) , 2 (2.4) (2.5) (2.6) f (xi1 ) + f (xi ) x, 2 i=1n

(2.3)

(2.7)

39 com Se = f (x0 ) + f (xn ), Si =n1 i=1

2.2 A regra do trapzio

(2.8) (2.9)

f (xi ).

Esta seo vale um mdulo de Python, que ns vamos denominar numint. Uma implementao razoavelmente eciente da regra do trapzio mostrada nas primeiras 21 linhas de numint.py, na listagen 2.5 Listagem 2.5: numint.py Integrao numrica, regra do trapzio1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n from __future__ import division def trapezio (n ,a ,b , f ): trapezio (n ,a ,b , f ): integra f entre a e b com n trapzios deltax = (b - a )/ n Se = f ( a ) + f ( b ) Si = 0.0 for k in range (1 , n ): xk = a + k * deltax Si += f ( xk ) I = Se + 2* Si I *= deltax I /= 2 return I # define Se # inicializa Si # calcula Si # clculo de I

O programa quadraver1.py calcula a integral de f (x) com 8 trapzios (listagem 2.6). Listagem 2.6: quadraver1.py Integrao numrica de f (x) com 8 trapzios1 2 3 4 5 6 7 8 9 10 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n from __future__ import division from numint import trapezio def f ( x ): return (( -1/48)* x **4 + (13/48)* x **3 + ( -17/12)* x **2 + (11/3)* x ) I3 = trapezio (8 ,1 ,5 , f ) print ( I3 = %8.4 f \ n % I3 )

A sada de quadraver1.py I3 14,6328. O erro est agora na 2a casa decimal, e o erro relativo = 0,0031, ou 0,3%. O problema com numint.trapezio que ns no temos uma idia do erro que estamos cometendo, porque, se estamos utilizando integrao numrica, porque no conhecemos o valor exato de I! Um primeiro remdio para este problema car recalculando a regra do trapzio com um nmero dobrado de trapzios, at que a diferena absoluta entre duas estimativas sucessivas que abaixo de um valor estipulado. Isto implementado, de forma muito ineciente, na prxima rotina

Matemtica Aplicada

40

do mdulo numint (listagem 2.7: note a continuao da numerao de linhas dentro do mesmo arquivo numint.py), denominada trapepsilonlento, e mostrada na listagem 2.7. Listagem 2.7: numint.py Integrao numrica ineciente, com erro absoluto prestabelecido20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 def t r a p e p s i l o n l e nt o ( epsilon ,a ,b , f ): t r a p e p s i l o n l en t o ( epsilon ,a ,b , f ): calcula a integral de f entre a e b com erro absoluto epsilon , de forma ineficiente eps = 2* epsilon n = 1 Iv = trapezio (1 ,a ,b , f ) while eps > epsilon : n *= 2 In = trapezio (n ,a ,b , f ) eps = abs ( In - Iv ) Iv = In return ( In , eps ) # # # # # # # # estabelece um erro inicial grande um nico trapzio primeira estimativa , " velha " loop dobra o nmero de trapzios estimativa " nova " , recalculada do zero calcula o erro absoluto atualiza a estimativa " velha "

O programa quadraver2.py calcula a integral de f (x) com erro absoluto estipulado menor que 0,0001, e imprime a estimativa da integral, o erro absoluto e o erro relativo (em relao ao valor exato conhecido) encontrados: I4 = 14,67777, = 0,00003, = 0,00000075. Com 4 casas decimais, este um resultado exato! Listagem 2.8: quadraver2.py Integrao numrica ineciente de f (x) com 0,00011 2 3 4 5 6 7 8 9 10 11 12 13 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n from __future__ import division from numint import t r a p e p s i l o n l en t o def f ( x ): return (( -1/48)* x **4 + (13/48)* x **3 + ( -17/12)* x **2 + (11/3)* x ) ( I4 , eps ) = t r a p e ps i l o n l e n t o (0.0001 ,1 ,5 , f ) print ( I4 = %8.5 f eps = %8.5 f % ( I4 , eps )) II = 1 4 . 6 7 7 7 7 7 7 7 7 7 7 7 7 7 7 delta = ( I4 - II )/ II print ( delta = %10.8 f % delta )

=

O problema com trapepsilonlento que todos os pontos que j haviam sido calculados por trapezio so recalculados em cada iterao (verique). Nosso ltimo esforo, a rotina trapepsilon em numint.py, corrige este problema, reaproveitando todos os clculos. Volte um pouco gura 2.2: nela, ns vemos a integrao numrica de f (x) com n = 4 e depois com n = 8 trapzios. Repare que Si para n = 4 parte do valor de Si para n = 8. De fato, para n = 4, Si = f (x2 ) + f (x4 ) + f (x6 ) (note que os ndices j esto denidos para o caso n = 8). Esta soma j foi calculada na integral com 4 trapzios, e no precisa ser recalculada. O que ns precisamos fazer agora somar f (x1 ) + f (x3 ) + f (x5 ) + f (x7 ) ao Si que j tnhamos. Se permanece o mesmo, e depois basta aplicar (2.7). Isto feito na ltima rotina de numint, denominada trapepsilon, na listagem 2.9.

41

2.2 A regra do trapzio

Listagem 2.9: numint.py Integrao numrica eciente, com erro absoluto prestabelecido35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 def trapepsilon ( epsilon ,a ,b , f ): trapepsilon ( epsilon ,a ,b , f ): calcula a integral de f entre a e b com erro absoluto epsilon , de forma eficiente eps = 2* epsilon n = 1 Se = f ( a ) + f ( b ) deltax = (b - a )/ n dx2 = deltax /2 Siv = 0.0 Iv = Se * dx2 while eps > epsilon : Sin = 0.0 n *= 2 deltax /= 2 dx2 = deltax /2 for i in range (1 , n ,2): xi = a + i * deltax Sin += f ( xi ) Sin = Sin + Siv In = ( Se + 2* Sin )* dx2 eps = abs ( In - Iv ) Siv = Sin Iv = In return ( In , eps ) # # # # # # # # # # # # # # # # # # # # estabelece um erro inicial grande n o nmero de trapzios Se no muda primeiro deltax primeiro deltax /2 Si " velho " I " velho " executa o loop pelo menos uma vez Si " novo " dobra o nmero de trapzios divide deltax por dois idem para dx2 apenas os mpares ... pula os ptos j calculados ! soma sobre os novos ptos internos aproveita todos os ptos j calculados I " novo " calcula o erro absoluto atualiza Siv atualiza Iv

O programa quadraver3.py (listagem 2.10) calcula a integral de f (x) com erro absoluto estipulado menor que 0,000001, e imprime a estimativa da integral, o erro absoluto e o erro relativo (em relao ao valor exato conhecido) encontrados: I5 = 14.6777776, = 0.0000005, = 0.00000001. Note que I5 exato at a 6a casa decimal, conforme estipulado. Listagem 2.10: quadraver3.py Integrao numrica eciente de f (x) com 0,0000011 2 3 4 5 6 7 8 9 10 11 12 13 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n from __future__ import division from numint import trapepsilon def f ( x ): return (( -1/48)* x **4 + (13/48)* x **3 + ( -17/12)* x **2 + (11/3)* x ) ( I5 , eps ) = trapepsilon (0.000001 ,1 ,5 , f ) print ( I5 = %12.7 f eps = %12.7 f % ( I5 , eps )) II = 1 4 . 6 7 7 7 7 7 7 7 7 7 7 7 7 7 7 delta = ( I5 - II )/ II print ( delta = %12.8 f % delta )

=

Exerccios Propostos 2.4 Usando Python e numint.trapezio, aproxime I = com 10 trapzios. 0

sen(x) dx pela regra do trapzio 0

2.5 Usando Python e numint.trapepsilon, aproxime I = trapzio com preciso absoluta menor ou igual a 1 105 .

sen(x) dx pela regra do

Matemtica Aplicada2.6 Dada f (x) denida no intervalo [x0 , x0 + 2h], deduza a regra de Simpson:x0 +2h

42

f (x) dx x0

h [f0 + 4f1 + f2 ] , 3

com fn = f (xn ), xn = x0 + nh, interpolando a funo g(x) = ax2 + bx + c atravs dos pontos (x0 , f0 ), (x1 , f1 ) e (x2 , f2 ) e calculando sua integral. 2.7 Dobrando o nmero de pontos de 2 para 4, obtenha a regra de Simpson para 5 pontos,x0 +4h

f (x) dx x0

h [f0 + 4f1 + 2f2 + 4f3 + f4 ] ; 3

generalize:bx0 +2nh

f (x) dx ax0

h [f0 + 4f1 + 2f2 + . . . + 2f2n2 + 4f2n1 + f2n ] . 3

2.8 Estenda numint com uma rotina simpson para calcular a regra de Simpson com 2n intervalos, e uma rotina simpepsilon para calcular uma integral numrica pela regra de Simpson com preciso estipulada. Baseie-se em trapezio e trapepsilon.

2.3 Aproximao de integrais com sries: a funo erroIntegraes numricas tais como as mostradas na seo 2.2 podem ser muito custosas em termos do nmero de operaes de ponto utuante necessrias. Algumas vezes, possvel ser mais inteligente. Um exemplo disto foi o clculo de I2 com uma nica parbola quadrtica1 no m da seo 2.1, que foi capaz de baixar o erro relativo para pouco mais de 1%. Considere agora o clculo de uma integral particularmente importante, a funo erro, denida por x 2 2 eu du. (2.10) erf(x) u=0 A funo erro est denida em Gnuplot e em Maxima, mas no em Python < 2.7. Como obt-la? Uma maneira fora bruta utilizar trapepsilon com uma preciso razovel (digamos, 106 ), gerar um arquivo, e plotar o resultado para ver a cara da erf(x). O programa vererf.py (listagem 2.11) calcula a funo em um grande nmero de pontos; gera um arquivo de dados vererf.dat, e ento o programa vererf.plt (listagem 2.12) plota o resultado, mostrado na gura 2.3. Observe que vererf.plt utiliza a funo erf pr-denida em Gnuplot e plota uma comparao: a linha contnua a erf de Gnuplot, e os pontos so os valores obtidos por vererf.py. Existe uma maneira mais inteligente de calcular erf(x): ela se baseia em integrar 2 a srie de Taylor do integrando, eu . Maxima permite calcular os primeiros termos:(% i1 ) taylor ( exp ( - u ^2) , u ,0 ,10) ; (% o1 )/ T / 1 - u 2 4 6 8 10 u u u u + -- - -- + -- - --- + . . . 2 6 24 120

1

Sim! existem parbolas cbicas, qurticas, etc..

43

2.3 Aproximao de integrais com sries: a funo erro

Listagem 2.11: vererf.py Clculo da funo erro por integrao numrica1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n from __future__ import division from math import exp , sqrt , pi from numint import trapepsilon def f ( x ): return exp ( - x * x ) xmin = 0.0 xmax = 3.0 nx = 100 dx = ( xmax - xmin )/ nx erf = 0.0 fou = open ( vererf . dat , wt ) xl = xmin fou . write ( %8.6 f %8.6 f \ n % ( xl , erf )) for k in range ( nx ): xu = xl + dx (I , eps ) = trapepsilon (1.0 e -6 , xl , xu , f ) erf = erf + (2.0/ sqrt ( pi ))* I fou . write ( %8.6 f %8.6 f \ n % ( xu , erf )) xl = xu fou . close ()

# # # # # # # # # # # # # # #

de 0 a 3 em 100 passos de dx erf (0) = 0 arquivo de sada limite inferior a partir de xmin erf (0) = 0 loop novo limite superior integra mais uma fatia acumula erf imprime at aqui atualiza limite inferior fecha o arquivo de sada

Listagem 2.12: vererf.plt Plotagem da funo erro por integrao numrica versus a erf(x) pr-denida em Gnuplot1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 set encoding iso_8859_1 # ------------------------------------------------------------------------------# no mundo de linux # ------------------------------------------------------------------------------set terminal postscript eps monochrome Times - Roman 18 set output vererf . eps set format y %3.1 f set grid set xlabel x set ylabel erf ( x ) set xrange [0:3] set yrange [0:1] set xtics 0 ,0.5 set ytics 0 ,0.1 plot vererf . dat using 1:2 title trapepsilon with points pt 7 ps 0.75 ,\ erf ( x ) title Gnuplot with lines lt 1 lw 2 # -----------------------------------------------------------------------------# agora no mundo de Windows # -----------------------------------------------------------------------------set terminal emf monochrome Times New Roman 18 set output vererf . emf replot exit

Matemtica Aplicada1.0 0.9 0.8 0.7 0.6 erf(x) 0.5 0.4 0.3 0.2 0.1 0.0 0 0.5 1 1.5 x 2 2.5 3 trapepsilon Gnuplot

44

Figura 2.3: Funo erf(x) calculada por integrao numrica, com trapepsilon e = 1 106 , versus a erf pr-denida em Gnuplot.

A expanso em torno de u = 0, e feita at o 10o termo. No muito difcil reconhecer nos denominadores os fatoriais de 0, 1, 2, 3, 4 e 5, e as potncias dos dobros destes valores nos expoentes de u. Em geral,

eu =

2

(1)nn=0

u2n . n!

(2.11)

Portanto,

x 0

eu du = = = =

2

x

(1)n n! n=0

0 n=0

(1)n0

u2n du n! u2n du

x

(1)n x2n+1 n! 2n + 1 n=0

(1)nn=0

x2n+1 . (2n + 1)n!

(2.12)

Cuidado: no toda srie que permite a troca impune da ordem de integrao e da soma (innita) da srie. Em princpio, devemos procurar os teoremas relevantes que nos permitem esta troca de ordem. Admitindo que est tudo bem, entretanto, ns conseguimos a srie de Taylor de erf(x)!

451.0 0.9 0.8 0.7 0.6 erf(x) 0.5 0.4 0.3 0.2 0.1 0.0 0

2.3 Aproximao de integrais com sries: a funo erroerf_1 Gnuplot

0.5

1

1.5 x

2

2.5

3

Figura 2.4: Funo erf(x) calculada com srie de Taylor, com erf_1, versus a erf pr-denida em Gnuplot. Ainda falta encontrar uma maneira computacionalmente eciente de passar do termo n 1 para o termo n. Vamos a isto: x2n+1 An (2n + 1)n! Bn Cn = = =

(1)n x2(n1+1)+1

2(n 1 + 1) + 1) n(n 1)! (1)n1 x2(n1)+1 [x2 ] 2(n 1) + 1) + [2] [n](n 1)! An1 (x2 ) (Bn1 + 2)(Cn1 n) (2.13)

Em outras palavras, o numerador An de cada novo termo da srie o anterior vezes x2 . O denominador mais complicado; ele formado pelo produto Bn Cn , e as relaes completas de recurso so An = x2n+1 = An1 (x2 ), Bn = 2n + 1 = Bn1 + 2, Cn = n! = Cn1 n. Mais uma coisa: a srie de erf(x) s contm potncias mpares: portanto, erf(x) uma funo mpar, e vale a relao erf(x) = erf(x). Com isto, temos uma rotina em Python para calcular erf(x), chamada erf_1(x), no mdulo erfs, do arquivo erfs.py. As linhas correspondentes a erf_1(x) so mostradas na listagem 2.13 Agora, podemos escrever vererf1.py (listagem 2.14), gerar um arquivo de sada vererf1.dat e plotar com o programa vererf1.plt (listagem 2.15). Os resultados continuam idnticos erf(x) de Gnuplot.

Matemtica Aplicada

46

Listagem 2.13: Clculo de erf(x) com uma srie de Taylor.1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n from __future__ import division from math import sqrt , pi def erf_1 ( x ): epsilon_erf = 1.0 e -6 eps_erf = 2* epsilon_erf A = x B = 1 C = 1 n = 0 termo = A /( B * C ) s = termo while eps_erf > epsilon_erf : n += 1 A *= ( - x * x ) B += 2 C *= n termo = A /( B * C ) eps_erf = abs ( termo ) s += termo return (2/ sqrt ( pi ) * s , n )

# # # # # # # # # # # # # # # #

mesma preciso garante entrada no while primeiro A primeiro B primeiro C primeiro n primeiro termo da srie primeira soma da srie loop incrementa n novo A novo B novo C novo termo seu valor absoluto soma na srie

Listagem 2.14: vererf1.py Clculo da funo erro com srie de Taylor entre 0 e 3.1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n from __future__ import division from math import exp , sqrt , pi from erfs import erf_1 xmin = 0.0 xmax = 3.0 nx = 100 dx = ( xmax - xmin )/ nx fou = open ( vererf1 . dat , wt ) xl = xmin fou . write ( %8.6 f %8.6 f \ n % ( xl ,0.0)) for k in range ( nx ): xu = xl + dx ( erf , n ) = erf_1 ( xu ) fou . write ( %8.6 f %8.6 f \ n % ( xu , erf )) xl = xu fou . close ()

# # # # # # # # #

de 0 a 3 em 100 passos de dx arquivo de sada limite inferior a partir de xmin erf (0) = 0 loop novo limite superior

# imprime at aqui # atualiza limite inferior # fecha o arquivo de sada

47

2.3 Aproximao de integrais com sries: a funo erro

Listagem 2.15: vererf1.plt Plotagem da funo erro calculada com srie de Taylor versus a erf(x) pr-denida em Gnuplot1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 set encoding iso_8859_1 # ------------------------------------------------------------------------------# no mundo de linux # ------------------------------------------------------------------------------set terminal postscript eps monochrome Times - Roman 18 set output vererf1 . eps set format y %3.1 f set grid set xlabel x set ylabel erf ( x ) set xrange [0:3] set yrange [0:1] set xtics 0 ,0.5 set ytics 0 ,0.1 plot vererf1 . dat using 1:2 title erf_1 with points pt 7 ps 0.75 ,\ erf ( x ) title Gnuplot with lines lt 1 lw 2 # -----------------------------------------------------------------------------# agora no mundo de Windows # -----------------------------------------------------------------------------set terminal emf monochrome Times New Roman 18 set output vererf1 . emf replot exit

Embora formalmente correta, a rotina erf_1 ca mais lenta medida em que |x| cresce. Como j vimos nas guras 2.3 e 2.4, erf(x) 1 para |x| 3. Porm, uma olhada em vererf1.dat, nos d erf_1(3) = 0.999978 , que ainda est acima da preciso especicada de 0,000001. Precisamos descobrir o valor de x para o qual erf_1 produz 0,999999. Por tentativa e erro, encontramos erf_1(3.6) = 0.99999952358847277 , que igual a 1,0 para 6 casas decimais. Precisamos tambm saber com quantos termos este clculo foi feito, e este o motivo de erf_1 devolver tambm n. Para calcular erf(3,6) com erf_1, ns precisamos de n = 43 termos. A rotina erf na listagem 2.16, tambm no arquivo erfs.py, usa este fato para impedir que o nmero de termos continue crescendo. Verique, voc mesmo(a), que erf(100) retorna 1.0, mas que erf_1(100) d um erro de ponto utuante. preciso enfatizar que erf em erfs.py ainda no uma rotina prossional. O nmero de termos usado ainda potencialmente muito grande; consequentemente, h muitas operaes de ponto utuante envolvendo A, B e C, e muitos testes dentro do while. possvel obter frmulas que aproximam erf(x) com menos termos, uniformemente, no intervalo (digamos) de 0 a 3,6: veja, por exemplo, Abramowitz e Stegun (1972). Mesmo assim, erf uma opo vantajosa em relao integrao numrica. A lio desta seo a seguinte: em geral, com esforo analtico adicional, possvel obter uma mistura de mtodos analticos e numricos que costuma ser amplamente superior ao uso de um mtodo numrico puro (por exemplo a regra do trapzio) para a obteno de resultados semelhantes. Exemplos deste fato vo reaparecer nos captulos seguintes.

Matemtica Aplicada

48

Listagem 2.16: erfs.py Clculo de erf(x) com srie de Taylor, limitado a no mximo 43 termos25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 def erf ( x ): if x > 3.6: return 1.0 elif x < -3.6: return -1.0 epsilon_erf = 1.0 e -6 eps_erf = 2* epsilon_erf A = x B = 1 C = 1 n = 0 termo = A /( B * C ) s = termo while eps_erf > epsilon_erf : n += 1 A *= ( - x * x ) B += 2 C *= n termo = A /( B * C ) eps_erf = abs ( termo ) s += termo return 2/ sqrt ( pi ) * s # limita o nmero de termos a 43

# # # # # # # # # # # # # # # #

mesma preciso garante entrada no while primeiro A primeiro B primeiro C primeiro n primeiro termo da srie primeira soma da srie loop incrementa n novo A novo B novo C novo termo seu valor absoluto soma na srie

Exerccios Propostos 2.9 (Bender e Orszag (1978), seo 6.2, Exemplo 2) Obtenha uma srie parax

F (x)

0

t1/2 et dt;

compare o resultado obtido com integrao numrica.

3Soluo numrica de equaes diferenciais ordinrias3.1 Soluo numrica de equaes diferenciais ordinriasSoluo analtica de uma EDO com Maxima. Considere uma equao diferencial de 1a ordem simples, forada eternamente por um seno: dy y + = sen x. dx x Na listagem 3.1, ns resolvemos esta equao com Maxima. Listagem 3.1: resolve-eqdif Soluo de uma EDO com Maxima1 2 3 4 5 6 7 8 9 10 (% i2 ) (% o2 ) (% i3 ) (% o3 ) y dy - + -- = sin ( x ) x dx dy y -- + - = sin ( x ) dx x ode2 (% , y , x ) sin ( x ) - x cos ( x ) + % c y = ---------------------x

(3.1)

Maxima nos informa de que a soluo geral da forma y(x) = sen x x cos x + c . x

evidente que, em geral, nem a equao diferencial nem sua soluo "existem"em x = 0. Entretanto, para c = 0, sen x y(x) = cos x. x Agora, sen x = 1, x de modo que existe uma soluo para a equao partindo de x = 0 se ns impusermos a condio incial y(0) = 0. De fato:x0

lim

x0

lim

sen x cos x = 1 1 = 0. x 49

Matemtica Aplicada1.5

50

1.0

0.5

y(x)

0.0

0.5

1.0

1.5

0

10

20x

30

40

50

Figura 3.1: Soluo da equao (3.1). O resultado est mostrado na gura 3.1. Claramente, existe uma parte transiente da soluo, dada por sen x/x, que morre medida que x cresce, e existe uma parte peridica (mas no permanente!) da soluo, dada por cos x, que domina y(x) quando x se torna grande. Ns dizemos que cos x parte estacionria da soluo. 3.1.1 Soluo numrica; mtodo de Euler A coisa mais simples que pode ser pensada para resolver a equao diferencial em questo transformar a derivada em uma diferena nita: y y + = sen x x x Isto um comeo, mas no suciente. Na verdade, o que desejamos que o computador gere uma lista de xs (uniformemente espaados por x), e uma lista de ys correspondentes. Obviamente, como os ys devem aproximar a funo, no podemos esperar deles que sejam igualmente espaados! Desejamos ento: x0 , x1 , . . . , xn onde com os correspondentes y0 , y1 , . . . , yn . Como x ser xo, podemos escrever nossa equao de diferenas nitas da seguinte forma: yi+1 yi y + = sen(x), x x onde eu deixei, propositadamente, x . . . = sen(x) y ainda sem ndices. De fato: qual xi e qual yi usar aqui? A coisa mais simples, mas tambm a mais instvel, usar i: yi+1 yi yi + = sen(xi ). x xi xi = ix,

51

3.1 Soluo numrica de equaes diferenciais ordinrias

Listagem 3.2: fracasso.py Um programa com o mtodo de Euler que no funciona1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* from __future__ import u n i c o d e _ l i t e ra l s from __future__ import print _functio n # ---------------------------------------------------------# cria listas com condies iniciais # ---------------------------------------------------------x = [0.0] y = [0.0] dx = 0.01 n = int (50/0.01) from math import sin for i in range (0 , n ): # de 0 at ... n -1 !!!! xn = ( i +1)* dx yn = y [ i ] + ( sin ( x [ i ]) - y [ i ]/ x [ i ])* dx x . append ( xn ) y . append ( yn ) fou = open ( fracasso . out , wt ) for i in range ( n ): fou . write ( %12.6 f %12.6 f \ n % ( x [ i ] , y [ i ])) fou . close ()

Note que agora possvel explicitar yi+1 em funo de todos os outros valores em i: yi+1 = yi + sen(xi ) yi x. xi

Este um exemplo de um mtodo explcito: o novo valor da funo em xi+1 (yi+1 ) s depende de valores calculados no valor anterior, xi . Um olhar um pouco mais cuidadoso ser capaz de prever o desastre: na frmula acima, se partirmos de x0 = 0, teremos uma diviso por zero j no primeiro passo! Muitos no veriam isto, entretanto, e nosso primeiro programa para tentar resolver a equao diferencial numericamente se chamar fracasso.py, e est mostrado na listagem 3.2 O resultado, mostrado na listagem 3.3 o seguinte fracasso: Listagem 3.3: Sada de fracasso.pyTraceback ( most recent call last ): File "./ fracasso . py " , line 15 , in < module > yn = y [ i ] + ( sin ( x [ i ]) - y [ i ]/ x [ i ])* dx Z e r o D i v i s i o n E r r o r : float division by zero

Isto j era previsvel: quando i == 0 no loop, x[0] == 0 no denominador, e o programa para com uma diviso por zero. Para conseguir fazer o mtodo numrico funcionar, ns vamos precisar de mais anlise! De volta equao diferencial, note que para que exista uma soluo na vizinhana de x = 0 necessrio que o limite limx0 y(x)/x exista; ns devemos ter y(x) = L, x0 x lim onde L uma constante nita a determinar. Mas, de nossa condio inicial,x0

lim y(x) = 0 L = 0.

Matemtica Aplicada

52

Listagem 3.4: sucesso.py Um programa com o mtodo de Euler que funciona1 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 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* # ---------------------------------------------------------# sucesso : resolve a equao diferencial # dy / dx + y / x = sen ( x ) # pelo mtodo de Euler . Uso : # # ./ sucesso . py < arquivo de sada > # ---------------------------------------------------------from __future__ import print _functio n from __future__ import division from sys import argv dx = float ( argv [1]) # tamanho do passo x = [0.0 , dx ] # condies inciais em x y = [0.0 , 0.0] # condies iniciais em y n = int (50/ dx ) from math import sin , cos # soluo numrica for i in range (1 , n ): xn = ( i +1)* dx yn = y [ i ] + ( sin ( x [ i ]) - y [ i ]/ x [ i ])* dx x . append ( xn ) y . append ( yn ) erro = 0.0 # calcula o erro relativo mdio for i in range (1 , n +1): yana = sin ( x [ i ])/ x [ i ] - cos ( x [ i ]) erro += abs ( ( y [ i ] - yana )/ yana ) erro /= n ; print ( erro relativo mdio = , %10.5 f % erro ) fou = open ( argv [2] , wt ) for i in range (0 , n +1): # imprime o arquivo de sada fou . write ( %12.6 f %12.6 f \ n % ( x [ i ] , y [ i ]) ) fou . close ()

Vamos ento reescrever a equao de diferenas usando o limite: y1 = y0 + sen(x0 ) y0 x = 0. x0

=0, lim x0 0

Na prtica, isto signica que ns podemos comear o programa do ponto x1 = x, y1 = 0! Vamos ento reescrever o cdigo, que ns agora vamos chamar, claro, de sucesso.py, que pode ser visto na listagem 3.4 A sada de sucesso.py gera o arquivo sucesso.out, que ns utilizamos para plotar uma comparao entre a soluo analtica e a soluo numrica, mostrada na gura 3.2 Na verdade, o sucesso estrondoso: com x = 0,01, ns conseguimos produzir uma soluo numrica que visualmente indistinguvel da soluo analtica. Uma das coisas que o programa sucesso.py calculou foi o erro relativo mdion

i=1

yi y(xi ) y(xi )

(onde yi a soluo numrica, e y(xi ) a soluo exata no mesmo ponto xi ). Para x = 0,01, = 0,02619, ou seja: menos de 3%. O preo, entretanto, foi alto: ns precisamos de um x bem pequeno, e de 50/0,01 = 5000 pontos para gerar a soluo. Ser possvel gerar uma soluo to boa com, digamos, 100 pontos?

53

3.1 Soluo numrica de equaes diferenciais ordinrias

1.5

x = 0,01 soluo numrica soluo analtica

1.0

0.5

y(x)

0.0

-0.5

-1.0

-1.5

0

10

20 x

30

40

50

Figura 3.2: Comparao da soluo analtica da equao (3.1) com a sada de sucesso.py, para x = 0,01.

1.5

x = 0,5 soluo numrica soluo analtica

1.0

0.5

y(x)

0.0

-0.5

-1.0

-1.5

0

10

20 x

30

40

50

Figura 3.3: Comparao da soluo analtica da equao (3.1) com a sada de sucesso.py, para x = 0,5.

Matemtica Aplicada

54

A gura 3.3 mostra o resultado de rodar sucesso.py com x = 0,5, muito maior do que antes. O erro mdio relativo agora pulou para = 1,11774, nada menos do que 111%, e muito pior do que a gura acima faz parecer primeira vista! 3.1.2 Um mtodo implcito, com tratamento analtico Nosso desao desenvolver um mtodo numrico que melhore consideravelmente a soluo mesmo com um x grosseiro, da ordem de 0,5. Nossa abordagem ser propor um mtodo implcito: xi + xi+1 yi+1 yi yi + yi+1 + = sen . x xi + xi+1 2 Note que tanto o termo y/x quando sen x esto sendo agora avaliados no ponto mdio entre xi e xi+1 . Lembrando que x = xi+1 xi , yi+1 : xi + xi+1 yi + yi+1 yi+1 yi = sen + , x xi + xi+1 2 1 1 1 1 xi + xi+1 + + yi + = sen , x xi + xi+1 x xi + xi+1 2 2xi+1 2xi xi + xi+1 yi+1 yi = sen . x(xi+1 + xi ) x(xi+1 + xi ) 2

yi+1

Uma rpida rearrumao produz yi+1 xi+1 yi xi = yi+1 x(xi+1 + xi ) xi + xi+1 sen , 2 2 xi x(xi+1 + xi ) xi + xi+1 = yi + sen . xi+1 2xi+1 2

(3.2)

Repare que a condio inicial y(0) = 0 no produz nenhuma singularidade em (3.2) para i = 0 x0 = 0, y0 = 0. O programa que implementa este esquema o sucimp.py, mostrado na listagem 3.5. O resultado um sucesso mais estrondoso ainda, e pode ser visto na gura 3.4. Agora, o erro mdio relativo baixou para = 0,02072, que ainda menor do que o do mtodo de Euler com x = 0,01, ou seja: com um x 50 vezes menor! 3.1.3 Runge-Kutta O preo que ns pagamos pelo mtodo extremamente acurado implementado em sucimp.py foi trabalhar analiticamente a equao diferencial, at chegar verso dedicada (3.2). Isto bom! Porm, s vezes no h tempo ou no possvel melhorar o mtodo por meio de um pr-tratamento analtico. Nossa discusso agora nos levar a um mtodo excepcionalmente popular, denominado mtodo de Runge-Kutta. Dada a equao dy = f (x, y) dx ns podemos tentar estimar a derivada no meio do intervalo h (estamos mudando a notao: at agora, usvamos x, mas h uma forma usual de explicitar o passo

55

3.1 Soluo numrica de equaes diferenciais ordinrias

Listagem 3.5: sucimp.py Mtodo de Euler implcito1 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 # !/ usr / bin / python # -* - coding : iso -8859 -1 -* # ---------------------------------------------------------# sucimp : resolve a equao diferencial # dy / dx + y / x = sen ( x ) # usando um mtodo implcito sob medida # ---------------------------------------------------------from __future__ import print _functio n from __future__ import division dx = 0.5 # passo em x x = [0.0] # x inicial y = [0.0] # y inicial n = int (50/ dx ) # nmero de passos from math import sin , cos # loop na soluo numrica for i in range (0 , n ): xn = ( i +1)* dx xm = x [ i ] + dx /2.0 yn = y [ i ]* x [ i ]/ xn + ( dx * xm / xn )* sin (( x [ i ]+ xn )/2) x . append ( xn ) y . append ( yn ) erro = 0.0 for i in range (1 , n