Método de Euler e solução de equações diferenciais
usando Fortran
Prof. Emílio Graciliano Ferreira Mercuri
Departamento de Engenharia Ambiental - DEA,
Universidade Federal do Paraná - UFPR
[email protected]
19 de fevereiro de 2013
Resumo: O objetivo da aula é apresentar uma introdução aos métodos aproximados para resolver
equações diferenciais ordinárias.
Introdução
O roteiro da aula deverá seguir os seguintes tópicos:
• Entrega dos Trabalhos em dupla
– Relatório em formato de artigo.
– Apresentar a modelagem do problema da evolução da temperatura do Lago.
• Cada dupla deverá portar o seu computador contendo:
– O compilador de Fortran Force.
– Arquivo com Aulas de Fortran.
• Introdução ao Método de Euler
– Método para resolução de EDO
– Aplicação de um exemplo em linguagem Fortran.
1
O Método de Euler (Método da Reta Tangente)
Considere um problema de valor inicial (PVI) representado por uma equação diferencial de
primeira ordem:
dy
= f (y, t)
dt
(1)
y(t0 ) = y0
(2)
e condição inicial
Sendo a Equação 1 linear e as funções f e fy contínuas em algum retângulo no plano
ty contendo o ponto (t0 , y0 ), existe um teorema no cálculo que garante uma solução única
y = φ(t) do problema dado em algum intervalo em torno de t0 .
O Método de Euler, ou método da reta tangente é expresso pela equação [?]:
yn+1 = yn + f (tn , yn )(tn+1 − tn ) n = 0, 1, 2, ...
(3)
Se o tamanho do passo tiver valor uniforme h e se denotarmos f (tn , yn ) = fn , então a
Equação 3 fica na forma:
yn+1 = yn + fn h n = 0, 1, 2, ...
(4)
O método de Euler consiste em calcular, repetidamente, a Equação 3 ou 4, usando o
resultado de cada passo para executar o próximo passo. Assim, obtemos a sequência de
valores
y0 , y1 , y2 , ..., yn , ...
(5)
que aproximam o valor da solução φ(t) nos pontos
t0 , t1 , t2 , ..., tn , ...
(6)
Um programa de computador para o método de Euler tem a estrutura a seguir. As
instruções específicas serão escritas na linguagem de programação Fortran. O algoritmo deverá
seguir os passos mostrados na Tabela 1.
Passo
Passo
Passo
Passo
Passo
Passo
1.
2.
3.
4.
5.
6.
Passo 7.
Passo 8.
Tabela 1: Algoritmo do Método de Euler
Defina f (t, y)
Insira os valores iniciais t0 e y0
Insira o tamanho do passo h e o número de passos n
Escreva t0 e y0
Para j de 1 até n Calcule
k1 = f (t, y)
y = y + h ∗ k1
t=t+h
Escreva t e y
Fim
Exemplo de aplicação em Fortran
Considere o PVI, para a função y = y(t) e sendo y 0 =
dy
:
dt
y 0 =1 − t + 4y
y(0) =1
(7)
A Eq. 7 é uma Equação linear de primeira ordem cuja solução que satisfaz as condições
iniciais pode ser encontrada em etapas. Começando pela equação homogênea:
dy
− 4y = 0
dt
(8)
dy
= 4dt
y
(9)
ln y = 4t + C1
(10)
y = Ce4t
(11)
A solução da equação não homogênea pode ser escrita na forma:
y = Ce4t + A + Bt
(12)
y 0 = 4Ce4t + B
(13)
4Ce4t + B − 4(Ce4t + A + Bt) = 1 − t
(14)
(B − 4A) − 4Bt = 1 − t
(15)
1
4
(16)
3
16
(17)
derivando:
e substituindo na EDO original:
−4B = −1
1
− 4A = 1
4
B=
A=−
e a Equação fica:
y = Ce4t −
3
1
+ t
16 4
(18)
3
=1
16
(19)
Aplicando a condição inicial y(0) = 1:
y(0) = C −
3
16 + 3
19
=
=
16
16
16
Portanto, a solução da equação diferencial é:
C =1+
1
3
19
y = t−
+ e4t
4
16 16
(20)
(21)
Como conhecemos a solução da EDO, a princípio não é necessário resolver o problema
numericamente. Porém, nem todas as equações possuem solução analítica e a disponibilidade
da solução exata torna fácil determinar a precisão de qualquer procedimento numérico utilizado
nesse problema.
A solução da EDO apresenta bastante divergência numérica dependendo do método escolhido. Agora vamos procurar qual o método mais preciso para resolvê-la, começando pelo
método de Euler.
Algoritmo
A seguir é apresentado o código do algoritmo na linguagem Fortran.
Code 1: Código em Fortran para aplicação do Algoritmo
1
PROGRAM EULER
2
3
4
5
6
7
8
REAL t , f
REAL y , h
REAL k1
REAL exata , e r r o
INTEGER i , n
!REAL VER
9
10
11
12
13
14
WRITE( ∗ , ∗ ) ’METODO DE EULER’
WRITE( ∗ , ∗ ) ’O PROGRAMA RESOLVE UMA EQUACAO DIFERENCIAL ORDINARIA’
15
16
! open ( 1 0 , f i l e = "EULER. t x t " )
17
18
19
20
t = 0.0
y = 1.0
h = 0.001
21
22
f = 1 − t + 4∗ y
23
24
n = 200 ! N mero de i t e r a es
25
26
DO i = 1 , n , 1
27
28
29
30
31
32
k1 = 1 − t + 4∗ y
y = y + h ∗ k1
e x a t a = 0 . 2 5 ∗ t −3.0/16.0 + 1 9 . 0 ∗ exp ( 4 . 0 ∗ t ) / 1 6 . 0
e r r o = abs ( exata−y )
WRITE( ∗ , ∗ ) t , y , exata , e r r o
33
34
35
! w r i t e ( 1 0 , 9 9 9 ) t , y , exata , e r r o
!999
format ( 4 F32 . 1 6 )
36
37
t = t + h
38
39
END DO
40
41
42
!VER = SYSTEM( ’ Wgnuplot u2 . gnu ’ )
!VER = SYSTEM( ’ OutPut . png ’ )
43
44
END PROGRAM EULER