Integração Numérica: Equações Diferenciais Pedro Barahona DI/FCT/UNL Introdução aos Computadores e à Programação 2º Semestre 2007/2008

Integração Numérica: Equações Diferenciais

Pedro BarahonaDI/FCT/UNL

Introdução aos Computadores e à Programação2º Semestre 2007/2008

9 Abril 2008 Integração Numérica: Equações Diferenciais Ordinárias 2

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!).

• 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).

dxdx

xdFxF

)()(

a x x+dx

f(x)

b

df(x)

dx



• 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.

a x x+h

f(x)

b

h

dA(x)

dx

xdA

h

xAhxAxf

h

)()()(lim)(

0



• 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.

bx

ax

bx

ax

bx

ax

bx

ax

ba hxfdxxfdt

dt

xdAdAA )()(

)(

a

f(x)

bh



• 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

• 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).

42)25(22)( 225

2

25

2

52

x

x

x

xxxFA

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; endforendfunction

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; endforendfunction


Cálculo de

• Exemplo: Para determinar o valor de , basta calcular a área da função

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).

21 xy

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

1

1

y1y2


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

dado instante. Denotando a posição do corpo no instante t por x(t) temos

• 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

5)(

dt

tdx



• 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

• 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

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.

...!3

)0(

!2

)0()0()0()(

3

3

32

2

2

h

dt

fdh

dt

fdh

dt

dffhf

dtdt

dffdtf

)0()0()(

erro

dt



• 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+1), 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; T(1) = 0; for i = 2:n 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

dtdt

dftfdttf

)0()()(



• 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

– a aceleração é a variação da velocidade em ordem ao tempo

• 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.

dt

tvdta

)()(

dt

txdtv

)()(


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


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

ai = - g

vi = vi-1 + ai-1 dt

yi = yi-1 + vi-1 dt


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; A(i) = - g - k*V(i-1); 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

ai = - g - k vi-1

vi = vi-1 + ai-1 dt

yi = yi-1 + vi-1 dt


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

• Já an horizontal, passamos a ter, as seguintes componentes

ayi = - g - k vyi-1

vyi = vyi-1 + ayi-1 dt

yi = yi-1 + vyi-1 dt

axi = - k vxi-1

vxi = vxi-1 + axi-1 dt

xi = xi-1 + vi-1 dt


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


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.

• Tudo contabilizado podemos modelar o movimento de uma mola através da

equação

ou seja, fazendo c = k/m e d = a/m

avkxdt

xdmF

2

2

x0

dvcxdt

xd

2

2


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; 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; endwhile t = T(i); Z = zeros(1,i); plot(T,[A;X;V;Z]);endfunction

dvcxdt

xd

2

2


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.


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.

logo ...!3

)0(

!2

)0()0()0()(

3

3

32

2

2

h

dt

fdh

dt

fdh

dt

dffhf

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

!2

)0()0()0()(

2

2

2 h

dt

fdh

dt

dffhf

Documents

Integração Numérica: Equações Diferenciais Pedro Barahona DI/FCT/UNL Introdução aos Computadores e à Programação 2º Semestre 2007/2008