Integração Numérica:
Equações Diferenciais
Pedro Barahona
DI/FCT/UNL
Introdução aos Computadores e à Programação
1º Semestre 2007/2008
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
1
Integração Numérica
•
Como é sabido, o integral de uma função nem sempre tem uma solução analítica
(priimitiva). Mesmo quando ela existe, pode não ser muito fácil determiná-la (derivar é
fácil, primitivar não!).
F ( x)  
•
dF ( x)
dx
dx
Nestas condições, podemos obter soluções aproximadas, através da integração
numérica. Estas aproximações são muito úteis no cálculo de áreas, em que a área
A(x) é a primitiva da função f(x).
dx
df(x)
f(x)
a
26 Outubro 2007
x x+dx b
Integração Numérica: Equações Diferenciais Ordinárias
2
Integração Numérica
•
Com efeito, denotando por A(x) a área da função f(x) desde a até x, e A(x+h) a área
da função entre a e x+h, a diferença destas áreas é dada por
A(x+h) – A(x) ≈ f(x) * h
•
i.e.
f(x) ≈ ( A(x+h)-A(x) ) / h
Esta aproximação só é exacta no limite, o que mostra ser a área A de uma função f a
primitiva dessa função (pois a derivada da área é a função, ou seja f(x) = dA(x)/dx.
f ( x)  lim
h 0
A( x  h)  A( x) dA( x)

h
dx
h
f(x)
dA(x)
a
26 Outubro 2007
x x+h
b
Integração Numérica: Equações Diferenciais Ordinárias
3
Integração Numérica
•
Assim, a obtenção da primitiva de uma função “corresponde” à obtenção da sua área.
Mas esta pode ser obtida numéricamente, substituindo-se o integral entre dois pontos
pelo somatório de áreas entre os dois pontos limites a e b, com incrementos dx = h.
x b
x b
x b
dA( x)
b
Aa   dA  
dt   f ( x)  dx 
dt
xa
x a
x a
x b
 f ( x)  h
xa
f(x)
a
26 Outubro 2007
h
b
Integração Numérica: Equações Diferenciais Ordinárias
4
Integração Numérica
•
Exemplo: Consideremos a função f(x) = 4x. Como neste caso, a função é facilmente
primitivável (F(x) = 2x2+k), a sua área entre os pontos x= 1 e x = 5 é dada por
x 5
A  F ( x) x  2  2 x
5
2
•
2 x 5
x 2
 2(52  22 )  42
Mas este é o valor aproximado obtido através das funções abaixo, em que a área é
aproximada com n intervalos (a função area mais geral, pois recebe o nome da
função a integrar como parâmetro).
function a =area_f(a,b,n);
s = 0;
d = (b-a)/n
for i = 1:n
y = 4*(a+d*(i-1));
s = s + y * d;
endfor
endfunction
26 Outubro 2007
function a = area(f,a,b,n);
s = 0;
d = (b-a)/n
for i = 1:n
y = feval(f,a+d*(i-1))
s = s + y * d;
endfor
endfunction
Integração Numérica: Equações Diferenciais Ordinárias
5
Cálculo de p
•
Exemplo: Para determinar o valor de p, basta calcular a área da função
y  1 x2
O valor de pi vem calculado por excesso em p1 e por defeito em p2. Em p vem a
média (correspondente à soma de trapézios, e não rectângulos).
function [p, p1, p2] = pi_integral(n);
s = 0; s1 = 0; s2 = 0;
x = 0; dx = 1/n;
for i = 1:n
y1 = sqrt(1-x^2);
x = x + dx;
y2 = sqrt(1-x^2);
s1 = s1 + y1 * dx;
s2 = s2 + y2 * dx;
s = s + 0.5 *(y1+y2) * dx;
endfor
p = 4*s; p1 = 4* s1; p2 = 4*s2;
endfunction
26 Outubro 2007
y1
y2
1
Integração Numérica: Equações Diferenciais Ordinárias
1
6
Equações Diferenciais Ordinárias
•
Uma outra aplicação da integração numérica ocorre na resolução numérica de
equações diferenciais.
•
Neste tipo de equações, não só são impostas restrições ao valor de variáveis., mas
também à sua variação.
•
Exemplo: Um corpo está sujeito a um movimento linear, ao longo do eixo dos xx,
com velocidade constante e igual a 5 m/s. Qual a posição desse corpo, sabendo-se
que no instante inicial ele está na posição x0 = 5m.
•
Resposta: Como é sabido, a velocidade instantânea mede a variação da posição num
addo instante. Denotando a posição do corpo no instante t por x(t) temos
dx (t )
5
dt
•
Neste caso a solução é simples e pode ser obtida analiticamente. Conhecendo-se a
derivada da distância, esta é obtida por primitivação da velocidade, o que neste caso
pode se feito analiticamente, obtendo-se
x(t) = 5 t + x0
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
7
Equações Diferenciais Ordinárias
•
Para uma dada função f, se se conhecer a sua variação, df/dt, ao longo do tempo,
e um valor inicial dessa variável, f0 = f(0), pode aproximar-se o valor f(t) da
variável ao longo do tempo. Com efeito, a função pode ser expandida numa série
de McLaurin (ou de Taylor) da função
df (0)
d 2 f (0) h 2 d 3 f (0) h3
f (h)  f (0) 
h

 ...
2
3
dt
dt
2!
dt
3!
•
Se o valor h for pequeno, os termos em h2, h3 e seguintes podem ser
desprezados. Neste caso, h é geralmente referido por dt (uma pequena variação
de t) e o valor da função no ponto t = h = dt é aproximado por
df (0)
f (dt )  f (0) 
dt
dt
erro
dt
para o que basta conhecer o valor de f no instante 0, f(0), a sua variação, df(0)/dt,
e o instante t = h = dt, suficientemente próximo de 0.
•
O processo pode depois repetir-se para sucessivos avanços dt até se obter o
valor de t pretendido.
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
8
Equações Diferenciais Ordinárias
•
Exemplo: Sabendo-se que a velocidade (isto é, a variação da posição) de um
corpo ao longo do tempo é dada por dx/dt = (3-x) + 8/t, o valor aproximado da
posição do corpo, num ponto t, x(t), pode ser obtido pela função abaixo, em que
– n é o número de intervalos dt considerados
– xo é a posição inicial
– xf, kx e kt são constantes (no caso, com valores 3, 1 e 8, respectivamente)
function X = corpo(x0, xf, kx, kt, t, n);
dt = t/n;
X(1) = x0
df (0)
T(1) = 0;
f
(
t

dt
)

f
(
t
)

dt
for i = 2:n
dt
T(i) = T(i-1)+ dt;
dxdt = kx*(xf-X(i-1)) + kt / (T(i-1)+1);
X(i) = X(i-1) + dxdt * dt;
endfor
plot(T,X);
endfunction
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
9
Equações Diferenciais Ordinárias
•
Em Física é muito frequente defrinirem-se grandezas como função não de outras
grandezas mas também da sua variação. Por exemplo,
– a velocidade é a variação da posição em ordem ao tempo
v (t ) 
d x (t )
dt
– a aceleração é a variação da velocidade em ordem ao tempo
a (t ) 
d v (t )
dt
•
Tendo em conta que a força é proporcional à aceleração, f = ma, o movimento de
um corpo (isto é, a função de posição ao longo do tempo x(t) ) pode ser
determinado a partir das forças que actuam sobre um corpo.
•
Vamos pois considerar alguns casos particulares de movimentos, quer numa só
dimensão quer a duas dimensões.
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
10
Equações do Movimento
•
Em geral, pretende-se determinar o valor da posição, x, velocidade, v, e
aceleração a, de um corpo num dado instante t, a partir da sua situação inicial
(velocidade inicial v0 e posição inicial x0) e das forças que sobre ele actuam.
•
Para esse efeito, iremos utilizar um conjunto de instantes t0, t1, t2, ..., tn, em que a
sua diferença dt = ti+1 – ti é um valor suficientemente pequeno para que os erros
produzidos pelas aproximações não sejam muito significativos.
•
Para facilitar a notação, representaremos por xi, vi e ai, respectivamente a
posição, velocidade e aceleração do corpo no instante ti. Nestas condições
– A partir da velocidade vi-1, e da sua variação, i.e. da aceleração ai-1 = dvi-1/dt,
obtem-se o valor da velocidade no instante seguinte
vi = vi-1 + ai-1 dt
– A partir da posição xi-1, e da sua variação, i.e. da velocidade vi-1 = dxi-1/dt,
obterm-se o valor da posição no instante seguinte
xi = xi-1 + vi-1 dt
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
11
Queda de Graves
•
Comecemos por considerar a queda de corpos sujeitos unicamente à força da
gravidade, que à superfície da terra tem um valor constante de 9.8 ms-2.
•
O seu movimento pode ser obtido pela função abaixo, em que v0 é a velocidade
inicial (positiva: para cima e negativa: para baixo) e y0 a posição inicial (se
positiva, representa posições acima do “chão”).
•
O movimento termina, em geral, quando o corpo atinge o chão (y <= 0).
function t = queda(y0, v0, dt);
i = 1;
T(1) = 0;
g = 9.8;
Y(1) = y0; V(1) = v0; A(1) = -g;
while Y(i) > 0
i = i+1;
A(i) = - g;
V(i) = V(i-1) + A(i-1) * dt;
Y(i) = Y(i-1) + V(i-1) * dt;
T(i) = T(i-1) + dt;
endwhile
t = T(i); Z = zeros(1,i); plot(T,[Y;Z]);
endfunction
26 Outubro 2007
ai = - g
vi = vi-1 + ai-1 dt
yi = yi-1 + vi-1 dt
Integração Numérica: Equações Diferenciais Ordinárias
12
Queda de Graves com Resistência do Ar
•
A resistência do ar, pode ser modelada como uma força proporcional à velocidade
e de sinal contrário, isto é, fa = -k v, em que k é um coeficiente de atrito de valor
constante. Assim no instante ti, teremos ai = - g + k*vi-1.
•
A função anterior pode ser alterada para reflectir esta mudança, como se indica
abaixo.
function t = queda_atrito(y0, v0, k, dt);
i = 1;
T(1) = 0;
g = 9.8;
Y(1) = y0; V(1) = v0; A(1) = -g - k*V(1);
while Y(i) > 0
i = i+1;
ai = - g - k vi-1
A(i) =
- g - k*V(i-1);
vi = vi-1 + ai-1 dt
V(i) = V(i-1) + A(i-1) * dt;
Y(i) = Y(i-1) + V(i-1) * dt;
yi = yi-1 + vi-1 dt
T(i) = T(i-1) + dt;
endwhile
t = T(i); Z = zeros(1,i); plot(T,[Y;Z]);
endfunction
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
13
Queda de Graves em 2D
•
Caso exista um velocidade inicial na horizontal, deverá ser considerado o
movimento em duas dimensões, com a anterior componente y, vertical e uma
nova componente x, horizontal.
•
Agora deverão considerar-se as duas componentes separadas, quer para a força
(aceleração), quer para a velocidade, quer ainda para a posição.
•
Na vertical, a situação mantem-se como anteriormente
ayi = - g
- k vyi-1
vyi = vyi-1 + ayi-1 dt
yi = yi-1 + vyi-1 dt
•
Já an horizontal, passamos a ter, as seguintes componentes
axi =
- k vxi-1
vxi = vxi-1 + axi-1 dt
xi = xi-1 + vi-1 dt
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
14
Queda de Graves em 2D
•
Desta forma a queda de graves em 2D pode ser simulada pela seguinte função
(em que o gráfico agora representa a evolução do corpo no espaço):
function t = queda_atrito_2D(x0, y0, vx0, vy0, k, dt);
i = 1;
T(1) = 0;
g = 9.8;
X(1) = x0; Vx(1) = vx0; Ay(1) =
- k*Vx(1);
Y(1) = y0; Vy(1) = vy0; Ay(1) = -g - k*Vy(1);
while Y(i) > 0
i = i+1;
Ay(i) =
- g
- k*Vy(i-1);
Vy(i) = Vy(i-1) + Ay(i-1) * dt;
Y(i) = Y(i-1) + Vy(i-1) * dt;
Ax(i) =
- k*Vx(i-1);
Vx(i) = Vx(i-1) + Ax(i-1) * dt;
X(i) = X(i-1) + Vx(i-1) * dt;
T(i) = T(i-1) + dt;
endwhile
t = T(i); Z = zeros(1,i); plot(X,[Y;Z]);
endfunction
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
15
Movimento de uma Mola
•
Um segundo exemplo de integração de equações diferenciais é obtido pelo
movimento de uma mola.
•
Se se deslocar o extremo da mola de uma distância x, ele é impelido para a posição
de repouso por uma força proporcional à distância. Existe ainda um factor de atrito,
proporcional e oposto à velocidade.
0
•
x
Tudo contabilizado podemos modelar o movimento de uma mola através da
equação
2
d x
F  m 2  kx  av
dt
ou seja, fazendo c = k/m e d = a/m
26 Outubro 2007
d 2x
 cx  dv
2
dt
Integração Numérica: Equações Diferenciais Ordinárias
16
Movimento de uma Mola
•
Desta forma o movimento da mola pode ser simulado pela função abaixo.
•
De notar o critério de paragem (velocidade baixa junto à posição de repouso, e
não sendo atingida, uma limitação no número de iterações).
•
De notar ainda a “instabilidade” criada se se usar um dt muito “grande”,
ocasionando oscilações crescentes!
function t = mola_atrito(x0, v0, c, d, dt, N);
i = 1;
T(1) = 0;
X(1) = x0; V(1) = v0; A(1) = -c*X(1) - d*V(1);
while i < N
i = i+1;
T(i) = T(i-1) + dt;
2
A(i) = -c*X(i-1) - d*V(i-1);
V(i) = V(i-1) + A(i-1) * dt;
2
X(i) = X(i-1) + V(i-1) * dt;
endwhile
t = T(i); Z = zeros(1,i);
plot(T,[A;X;V;Z]);
endfunction
d x
 cx  dv
dt
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
17
Movimento de uma Mola
•
Em algumas funções, nomeadamente as periódicas, os pequenos erros
acumulados podem causar importantes desvios entre a função aproximada e a
função “real”.
•
Naturalmente, quanto menores forem os “dt” menores, em princípio serão os erros
nas simulações.
•
No entanto, este “remédio” solução levanta alguns problemas, nomeadamente:
1. As simulações tornam-se muito lentas, pois envolvem ciclos com maior
número de iterações e portanto um maior número de operações;
2. Os erros de arredondamento podem começar a ter uma grande influência no
resultado e eliminar a maior precisão obtida com menores incrementos do
tempo, dt, ou de outras variáveis de iteração.
•
Uma alternativa à diminuição do incremento da variável de iteração é a utilização
de melhores aproximações no cálculo da função.
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
18
Movimento de uma Mola
•
Em particular, sendo conhecida derivadas de ordem superior, podem estas ser
usadas na aproximação de uma função em séreis de Taylor, em particular as
segundas derivadas, que estão calculadas em alguns problemas anteriores.
df (0) d 2 f (0) h 2 d 3 f (0) h3
f (h)  f (0) 
h

 ...
dt
dt2 2!
dt3 3!
logo
df (0)
d 2 f (0) h 2
f (h)  f (0) 
h
dt
dt 2 2!
function t = mola_atrito_2(x0, v0, c, d, dt, N);
i = 1;
T(1) = 0;
X(1) = x0; V(1) = v0; A(1) = -c*X(1) - d*V(1);
while i < N
i = i+1;
T(i) = T(i-1) + dt;
A(i) = -c*X(i-1) - d*V(i-1);
V(i) = V(i-1) + A(i-1) * dt;
X(i) = X(i-1) + V(i-1) * dt + 1/2 * A2(i-1) * dt^2;
endwhile
t = T(i); Z = zeros(1,i);
plot(T,[A;X;V;Z]);
endfunction
26 Outubro 2007
Integração Numérica: Equações Diferenciais Ordinárias
19