Trabalho de Práticas Computacional de Mecânica

1

Trabalho de Práticas

Computacionais de Mecânica

Hugo Torres - 120322043 João Guerra - 120302025 Liliana Andrade - 120322039 Simão João - 120302050 Simão Sá – 120302048

Data: 28 de Novembro 2012

2

Índice 1.Problema 1 .............................................................................................................................................. 4

1.1.Introdução ......................................................................................................................................... 4

1.2.Determinação da solução exacta ...................................................................................................... 4

1.3.Aplicação do método de Euler .......................................................................................................... 5

1.4.Como determinar o erro absoluto? .................................................................................................. 6

1.5.Será que a dependência do erro é linear? ........................................................................................ 7

1.6.Será que a dependência do erro é quadrática? ................................................................................ 8

1.7.Conclusão .......................................................................................................................................... 9

1.8.Análise do código do Problema 1 ...................................................................................................10

2.Problema 2 ............................................................................................................................................12

2.1.Introdução .......................................................................................................................................12

2.2.Determinação das solução exactas .................................................................................................12

2.3.Comparação dos gráficos ................................................................................................................13

2.4.Conclusão ........................................................................................................................................14


3.Problema 3 ............................................................................................................................................17

3.1.Introdução .......................................................................................................................................17

3.2.Determinação da solução exacta ....................................................................................................17

3.3.Aplicação do Método de Runge-Kutta de segunda ordem .............................................................17

3.4.Como determinar o erro absoluto? ................................................................................................19

3.5.Comparação dos gráficos das aproximações ..................................................................................20

3.6.Conclusão ........................................................................................................................................20


4.Problema 4 ............................................................................................................................................22

4.1.Introdução .......................................................................................................................................22

4.2.Integração da equação diferencial ..................................................................................................23

4.3.Gráficos obtidos ..............................................................................................................................23

4.4.Conclusão ........................................................................................................................................24


3

5.Problema 5 ............................................................................................................................................26

5.1.Introdução .......................................................................................................................................26

5.2.Construção e análise dos gráficos ...................................................................................................26

5.3.Conclusão ........................................................................................................................................29

5.4.Análise do código do Problema 5……………………………………………………………………………………………….29

6.Problema 6 ............................................................................................................................................31

6.1.Introdução .......................................................................................................................................31

6.2.Comparação dos integradores ........................................................................................................32

6.3.Comparação e análise dos gráficos .................................................................................................33

6.4.Conclusão ........................................................................................................................................36


7.Problema 7 ............................................................................................................................................39

7.1.Introdução .......................................................................................................................................39

7.2.Como passar os argumentos para a lista de valores t? ..................................................................39

7.3.Testando o código ...........................................................................................................................40

7.4.Conclusão ........................................................................................................................................41


4

1. Problema 1

No problema de atrito de Stokes a solução exacta é conhecida. Estude a dependência no erro

de cálculo de v(t+h) em função de h num passo de Euler. A dependência é linear, quadrática?

1.1. Introdução

Este problema pede-nos para encontrar a relação entre o erro num passo do Método de Euler e

o tamanho deste passo quando este assume valores pequenos.

1.2. Determinação da solução exacta:

A força de atrito de stokes é uma força que depende da velocidade a que o corpo se desloca,

segundo a equação:

Integrando, obtemos:

∫

∫

|

|

Tomando :

5

1.3.Aplicação do Método de Euler:

Da definição de derivada, vem:

|

Para um h suficientemente pequeno:

|

Rearranjando os termos, obtemos um passo do Método de Euler:

|

|

Recordemos que, da equação diferencial, vem:

|

Logo, com , temos:

Para :

Este é o passo do Método de Euler para o instante inicial. Para obter os pontos seguintes, é

necessário repetir o passo , actualizando de cada vez que

queremos computar o próximo ponto. Como isto se trata de uma aproximação, é natural que

haja um erro subjacente, que é o que queremos estudar. Para efectuar no Python um passo do

Método de Euler, usamos o código da função eulerstep que nos foi fornecida nas aulas

práticas:

def eulerstep(func,x0,h):

return x0+func(x0)*h

6

1.4. Como determinar o erro absoluto?

Para encontrar o erro absoluto associado ao Método de Euler em função do tamanho do passo,

procedemos do seguinte modo:

1. Fixamos um ponto . Escolhemos o ponto inicial , pois já conhecemos o valor

da velocidade nesse ponto, .

2. Utilizando um passo do Método de Euler para um determinado h, determinamos o valor

aproximado da velocidade no ponto seguinte, .

3. Usando a solução exacta da equação que obtivemos, calculamos a velocidade nesse

mesmo ponto .

4. Achamos a diferença entre os dois valores , de modo a

encontrar o erro absoluto para esse novo ponto.

De modo a estudar a dependência do erro absoluto, repetimos este processo para vários valores

de h.

O gráfico obtido para o erro absoluto em função de h, através da função plot, é o seguinte:

erro=zeros((npoints,2),float)

i=0

for j in range(npoints):

erro[j,0]=i

erro[j,1]=abs(fsol(i)-eulerstep(fint,x0,i))

i=i+soma

i=0

for j in range(101):

erro[j,0]=i


i=i+soma

plot(erro[:,0],erro[:,1],'r')

xlabel('h')

ylabel('erro')

title('Estudo do erro num passo do Metodo de Euler')

legend(('Exacto'))

show()

7

1.5.Será que a dependência do erro é linear?

Tentamos agora encontrar a recta que melhor se ajusta aos dados. Isto é conseguido através de

uma regressão linear. No python, utilizamos a função polyfit (dentro de uma função

adjust que diz qual o grau do polinómio pretendido) para determinar os coeficientes dessa

recta. De seguida, criamos uma matriz, matriz, para armazenar os valores da recta para cada

valor do passo do Método de Euler, de modo a poder compará-los mais tarde.

Novamente, através da função plot, podemos obter este gráfico e a sobreposição dele com o

do erro absoluto exacto:

Facilmente se conclui, pela análise dos dois gráficos, que a dependência não é linear.

def adjust(ordem):

matriz=zeros((npoints,2),float)

coeford=polyfit(erro[:,0],erro[:,1],ordem)

i=0

def eq(x):

eq=0

for m in range(ordem+1):

eq=eq+coeford[m]*(x**(ordem-m))

m=m+1

return eq


matriz[j,0]=i

matriz[j,1]=eq(i)

i=i+next

return matriz

lin=adjust(1)

plot(erro[:,0],erro[:,1],'r',lin[:,0],lin[:,1],'g')

xlabel('h')

ylabel('erro')


legend(('Exacto','1 Grau'))

show()

8

1.6.Será que a dependência do erro é quadrática?

Como tentámos ajustar uma recta ao gráfico e falhámos, tentamos desta vez ajustar uma curva.

Utilizando o mesmo processo para armazenar os dados, recorremos de novo à função polyfit e

determinamos os coeficientes do polinómio de grau 2 que melhor se ajusta aos dados.

A função plot devolve-nos os seguintes gráficos:

Daqui podemos concluir que a dependência do erro num passo do Método de Euler é

quadrática, porque ajusta-se muito bem ao gráfico de um polinómio do segundo grau.

def adjust(ordem):



i=0

def eq(x):

eq=0



m=m+1

return eq


matriz[j,0]=i

matriz[j,1]=eq(i)

i=i+next

return matriz

quad=adjust(2)

plot(erro[:,0],erro[:,1],'r',quad[:,0],quad[:,1],'b--')

xlabel('h')

ylabel('erro')


legend(('Exacto','2 Grau'))

show()

9

1.7.Conclusão:

Utilizando o Método de Euler na equação diferencial que traduz o atrito de Stokes, sabemos

que há um erro associado à escolha do tamanho do passo. Este erro pode ser minimizado

tomando um passo muito pequeno, mas está sempre presente, pois este método trata-se de uma

aproximação. Assim, decidimos ver de que forma é que este erro varia com a escolha do passo

para termos uma melhor noção de como este método funciona. Ajustando rectas e curvas ao

gráfico do erro absoluto, chegámos à conclusão que a dependência do erro é quadrática, pois

ajusta-se ao gráfico de uma função polinomial do 2º grau.

10

1.8. Análise do código do Problema 1:

Começamos por importar os módulos do Python que nos interessam: scipy, pylab e

numpy, e definimos as variáveis que serão relevantes para o problema.

Em seguida, definimos a função eulerstep, sendo func a função a integrar, x0 o ponto da

curva solução e h o tamanho do passo do Método de Euler.

fsol é a solução exacta para a equação diferencial do atrito de stokes.

fint é a função a integrar numericamente. Vem da equação diferencial que

|

.

Definimos a matriz erro de duas colunas para armazenar os valores a por no gráfico.

Preenchemos a primeira coluna com o valor de h e segunda coluna com a diferença entre o

valor exacto e o valor aproximado obtido pelo Método de Euler.

from scipy import *

from pylab import *

from numpy import *

t0=0

x0=1

hmax=0.4

npoints=100

next=hmax/npoints



def fsol(x):

return exp(-x)

def fint(x):

return -x


i=0


erro[j,0]=i


i=i+next

i=0


erro[j,0]=i


i=i+soma

11

Depois, de modo a não ter de repetir o processo para cada grau de aproximação, criamos a

função adjust que determina o grau da função polyfit e nos devolve uma matriz com os

valores da recta/curva correspondente, para cada valor de h.

Assim, depois de criada esta função, torna-se fácil obter as matrizes lin e quad para uma

regressão linear e para uma regressão quadrática, respectivamente.

Finalmente, implementamos a função que cria os gráficos.

def adjust(ordem):



i=0

def eq(x):

eq=0



m=m+1

return eq


matriz[j,0]=i

matriz[j,1]=eq(i)

i=i+next

return matriz

lin=adjust(1)

quad=adjust(2)

plot(erro[:,0],erro[:,1],'r',lin[:,0],lin[:,1],'g',quad[:,0],

quad[:,1],'b--')

xlabel('h')

ylabel('erro')


legend(('Exacto','1 Grau','2 Grau'))

show()

12

2. Problema 2

Considere o caso de atrito de Newton,

Compare os amortecimentos de Stokes e Newton para valores iniciais de v e

iguais.

2.1.Introdução:

Este problema pede-nos para comparar dois regimes diferentes: o de atrito de Newton e o de

atrito de Stokes. No caso do atrito de Stokes, a força varia linearmente com a velocidade, mas

no caso do atrito de Newton, a força varia com o quadrado da velocidade, por isso é de esperar

que as funções se comportem de modo diferente.

2.2.Determinação das soluções exactas:

De modo a comparar analiticamente as duas funções, e ter uma referência quando

implementarmos o Método de Euler, faz sentido encontrar a solução exacta para cada uma

delas.

A solução para a força de atrito de Stokes já é conhecida do Problema 1:

E a solução para a força de atrito de Newton pode ser facilmente derivada:

∫

∫

|

|

13

Tomando as mesmas condições iniciais do problema anterior, .

Como a função do atrito de Stokes varia com uma exponencial e a do atrito de Newton com um

inverso do tempo, será de esperar que a primeira função decresça mais rapidamente do que a

segunda.

2.3.Comparação dos gráficos:

Para comparar os dois amortecimentos para valores de velocidade e aceleração iguais, tomamos

os coeficientes como parâmetros e calculamos a velocidade inicial em função destes.

Logo,

Podemos retirar o caso em que pois esse é o caso trivial em que os corpos não se

movem. Assim:

Tomando

e , podemos finalmente comparar os gráficos através da função

plot.

14

Vemos claramente que sob as mesmas condições iniciais e para coeficientes iguais, a função do

atrito de Stokes decresce mais rapidamente do que o do atrito de Newton, desde o instante

inicial, tal como previsto pela análise das soluções das equações diferenciais. Mas será que este

resultado se mantém para diferentes coeficientes?

Se olharmos simplesmente para o quociente entre as duas velocidades, nada nos garante qual

dos amortecimentos é mais forte no início, mas se tomarmos o limite quando , é fácil

ver que a exponencial vai dominar, independentemente do valor dos coeficientes. No entanto,

se os coeficientes forem iguais, fazendo uma expansão de Taylor para a exponencial à volta do

ponto , com

, ficamos com:

O que prova que para coeficientes e massas iguais, uma partícula sujeita ao amortecimento de

Stokes terá sempre velocidade menor do que uma partícula sujeita ao amortecimento de

Newton.

2.4.Conclusão:

O atrito de Stokes e o atrito de Newton dependem ambos da velocidade a que o corpo em

questão se move. No entanto, como o atrito de Stokes varia linearmente com a velocidade, vai

decrescer mais rapidamente do que o atrito de Newton, que varia com o quadrado da

velocidade. Quer isto dizer que uma partícula sujeita somente ao atrito de Stokes, parará mais

rapidamente do que uma partícula sujeita ao atrito de Newton.

plot(mstokes[:,0],mstokes[:,1],'r',mnewton[:,0],mnewton[

:,1],'g')

xlabel('tempo')

ylabel('velocidade')

title('Amortecimento de Stokes e de Newton')

legend(('Stokes','Newton',))

show()

15

2.5.Análise do código do Problema 2:

Começamos por importar os módulos que serão relevantes para o problema: scipy e pylab e

definimos o passo do Método de Euler em função da quantidade de pontos que queremos.

Em seguida, definimos a função eulerstep, sendo func a função a integrar, x0 o ponto da

curva solução e h o tamanho do passo do Método de Euler. A função euler preenche as

tabelas de dados com os pontos encontrados pela função eulerstep. t0 é o instante em

questão.

Definimos as funções stokes e newton que serão chamadas para a função eulerstep

pela função euler de modo a preencher as tabelas mstokes e mnewton, respectivamente.

From scipy import *

from pylab import *

npoints=100

h=3.0/npoints



def euler(func):

t0=0

x0=1

inserir=zeros((npoints+1,2),float)

for j in range(npoints+1):

inserir[j,0]=t0

inserir[j,1]=x0

x0=eulerstep(func,x0,h)

t0=t0+h

j=j+1

return inserir

def stokes(v):

return -v

def newton(v):

return -v**2

mstokes=euler(stokes)

mnewton=euler(newton)

16

Por fim, implementamos a função plot que nos devolve o gráfico.

plot(mstokes[:,0],mstokes[:,1],'r',mnewton[:,0],mnewton[

:,1],'g')

xlabel('tempo')


title('Amortecimento de Stokes e de Newton')

legend(('Stokes','Newton',))

show()

17

3. Problema 3

Resolva a mesma questão do Problema 1 na página 6 para um passo de Runge-Kutta de

segunda ordem.

3.1.Introdução

À semelhança do Problema 1, é-nos pedido para encontrar a dependência do erro em função de h

quando este assume valores pequenos, mas desta vez, usamos um método um pouco diferente, o Método

de Runge-Kutta de segunda ordem. Tal como o Método de Euler, obtemos somente uma aproximação

da função-solução exacta, portanto obteremos um erro, e é esse erro que queremos estudar.

3.2.Determinação da solução exacta:

A solução exacta para a equação diferencial do atrito de Stokes já foi determinada no Problema

1:

3.3.Aplicação do Método de Runge-Kutta de segunda ordem:

O Método de Runge-Kutta de ordem 2 assume que o declive da

recta que une os pontos e é aproximadamente igual

ao declive da recta tangente ao gráfico da função no ponto entre

e , o ponto . Assim, vem que:

|

|

|

Como

:

(

)

Como não sabemos o valor da função no ponto

, recorremos a uma expansão de Taylor

sobre o ponto :

18

Como este método é de segunda ordem, só nos interessam os dois primeiros termos:

(

)

(

)

(

)

Juntando as duas equações, obtemos:

(

)

Este é um passo do Método de Runge-Kutta de segunda ordem para a equação diferencial que

traduz o atrito de Stokes:

Agora, é-nos possível computar cada um dos pontos seguintes, dentro do erro desejado, com

um passo Método de Runge-Kutta de ordem 2. Mais uma vez, como se trata apenas de uma

aproximação, há um erro associado, que é o que queremos estudar. Para implementar esta

função no Python, usaremos a função rk2step que foi fornecida nas aulas práticas:

def rk2step(func, x0,h):

k0= func(x0)*h

return x0+func(x0+0.5*k0)*h

19

3.4.Como determinar o erro absoluto?

Tal como para o Método de Euler no Problema 1, para encontrar o erro absoluto associado ao

Método de Runge-Kutta de ordem 2 (RK2) em função do tamanho do passo, procedemos do

seguinte modo:

5. Fixamos um ponto . Escolhemos o ponto inicial , pois já conhecemos o valor

da velocidade nesse ponto, .

6. Utilizando um passo de RK2 para um determinado h, determinamos o valor aproximado

da velocidade no ponto seguinte, .

7. Usando a solução exacta da equação que obtivemos, calculamos a velocidade nesse

mesmo ponto .

8. Achamos a diferença entre os dois valores , de modo a

encontrar o erro absoluto para esse novo ponto.

De modo a estudar a dependência do erro absoluto, repetimos este processo para vários valores

de h.

O gráfico do erro em função de h é o seguinte:


i=0


erro[j,0]=i

erro[j,1]=abs(fsol(i)-rk2step(fint,x0,i))

i=i+soma

i=0


erro[j,0]=i


i=i+soma

plot(erro[:,0],erro[:,1],'b')

xlabel('h')

ylabel('erro')

title(‘Estudo do erro num passo do Metodo de Runge-Kutta

de segunda ordem’)

legend(('Exacto'))

show()

20

3.5.Comparação dos gráficos das aproximações:

Se sobrepusermos os gráficos, vemos claramente que a dependência do erro é cúbica, pois é o

gráfico da função polinomial cúbica que melhor se ajusta ao erro absoluto em função de h.

3.6.Conclusão:

Utilizando o Método de Runge-Kutta de segunda ordem na equação diferencial que traduz o

atrito de Stokes, podemos concluir que, embora ainda haja um erro subjacente ao método, este

erro é menor do que o erro obtido através do Método de Euler (basta olhar para o eixo das

ordenadas para ver que os valores do erro são muito mais pequenos). Isto permite-nos obter

uma aproximação numérica muito precisa e que não requer tanto tempo de processamento do

programa como o Método de Euler. Para além disso, verificamos que a dependência do erro

devido à escolha do valor de h é cubica. Isto deve-se ao facto de termos usado a expansão de

Taylor somente até ao termo do segundo grau.

21

3.7.Análise o código do Problema 3:

O código para o Problema 3 é essencialmente igual ao do Problema 1. As únicas diferenças

consistem no seguinte:

Utilizamos um integrador diferente. Em vez de eulerstep, usamos rk2step:

Depois de definida a função adjust, criamos mais uma matriz para além de lin e quad, a

matriz cube, que corresponde a uma regressão cúbica:

Finalmente, adicionamos um termo à função plot, correspondente a essa regressão cúbica:

def rk2step(func, x0,h):

k0=func(x0)*h


lin=adjust(1)

quad=adjust(2)

cube=adjust(3)

c

plot(erro[:,0],erro[:,1],'r',lin[:,0],lin[:,1],'g--',

quad[:,0],quad[:,1],'y--',cube[:,0],cube[:,1],'b--')

xlabel('h')

ylabel('erro')

title(Estudo do erro num passo do Metodo de Runge-Kutta

de ordem 2')

legend(('Exacto','1 Grau','2 Grau','3 Grau'))

show()

22

4. Exercício 4

Estude o caso de queda livre de um corpo na atmosfera (sentido positivo z, ascendente) pelo

método de rk2, determinando a função . A força segundo z é:

Considere queda na vertical, com velocidade inicial e produza gráficos

comparativos para diferentes valores de . Não se esqueça de reduzir dimensionalmente as

equações.

4.1.Introdução

Neste problema, é-nos pedido para estudar o caso de um corpo em queda livre sujeito ao

amortecimento de Newton através da construção de gráficos para vários valores do coeficiente

.

A equação que nos dá a força resultante no corpo é:

Como o movimento é sempre descendente, a velocidade terá sinal negativo, logo podemos

simplificar a equação:

Para obtermos a aceleração, basta dividir a equação por m:

Como o que nós queremos estudar é a dependência entre velocidade e a constante de Newton,

temos que efectuar uma substituição que tire as outras constantes da equação:

√

Deste modo temos:

√

Logo, teremos também que tomar um em função de t:

√

Obtendo-se assim o resultado final:

23

Com esta equação, podemos começar a estudar a dependência entre a velocidade e .

Adicionalmente, como sabemos que , então .

4.2. Integração da equação diferencial

Como já é explicado no problema será utilizado o Método de Runge-Kutta de ordem 2 para fazer a

aproximação da função. Neste trabalho irão ser utilizados 100 pontos, um h de 0,1 e 4 valores distintos

de : 0.33,1,4 e 8.

Com vista a agrupar as matrizes que tinham os dados sobre os tempos e as respetivas velocidades para

cada constante, criámos uma lista chamada matrizes onde foram colocadas as 4 matrizes

correspondentes aos 4 valores de .

4.3. Gráficos obtidos:

matrizes=[]

for gamma in[0.33,1,4,8]:

traj=zeros((100,2),float)

for i in range(100):

traj[i,0]=t0

traj[i,1]=x0

x0=rk2(az,x0,h)

t0+=h

matrizes.append(traj)

t0=0

x0=0

plot(matrizes[0][:,0],matrizes[0][:,1],'r',matrizes[1][:

,0],matrizes[1][:,1],'g',matrizes[2][:,0],matrizes[2][:,

1],'b',matrizes[3][:,0],matrizes[3][:,1],'y')

xlabel('tempo')


title('Variacao da velocidade de um corpo em queda livre

com resistencia do ar')

legend(('gama=0.33','gama=1','gama=4','gama=8'))

show()

24

4.4. Conclusão

Através da análise do gráfico obtido, concluímos que a velocidade terminal do corpo é tanto

maior quanto menor for . Isto deve-se ao facto de ser proporcional à resistência do ar

(como se pode ver na fórmula que nos é dada no problema), o que significa que quanto menor

for a resistência do ar, maior será a velocidade terminal atingida

Além disso, também se verifica que quanto maior for menos acelerado é o movimento. Isto

deve-se ao facto descrito no parágrafo anterior, porque quanto maior for , maior é a força

ascendente sobre o corpo. Assim, como a aceleração exercida pela resistência do ar no corpo é

maior, então a aceleração final será menor.

4.5. Análise do código do Problema 4:

Começamos por importar os módulos que nos interessam e definir as variáveis que serão

relevantes para o problema: o tempo inicial t0, a posição inicial x0 e o tamanho do passo h.

Em seguida, definimos a função az que nos dá a função a integrar:

rk2 é a função que representa o passo do Método de Runge-Kutta de segunda ordem.

from scipy import *

from pylab import *

from numpy import *

from matplotlib import *

t0=0

x0=0

h=0.1

def az(vz):

return gamma*(vz**2)-1

def rk2(func,x0,h):

k= func(x0)*h

return x0+func(x0+k/2)*h

25

Depois, criamos um ciclo para preencher as tabelas com os valores das funções para cada .

Por fim, usamos a função plot para desenhar o gráfico:

matrizes=[ ]

for gamma in[0.33,1,4,8]:

traj=zeros((100,2),float)

for i in range(100):

traj[i,0]=t0

traj[i,1]=x0

x0=rk2(az,x0,h)

t0+=h

matrizes.append(traj)

t0=0

x0=0

plot(matrizes[0][:,0],matrizes[0][:,1],'r',matrizes[1][:

,0],matrizes[1][:,1],'g',matrizes[2][:,0],matrizes[2][:,

1],'b',matrizes[3][:,0],matrizes[3][:,1],'y')

xlabel('tempo')


title('Variacao da velocidade de um corpo em queda livre

com resistencia do ar')

legend(('gama=0.33','gama=1','gama=4','gama=8'))

show()

26

5. Problema 5

O problema 2 na página II ficou um pouco incompleto por não termos determinado em

conjunto com . Corrija este facto usando um sistema de equações de primeira ordem. Estude de

novo o caso de queda livre de um corpo na atmosfera (sentido positivo z, ascendente) pelo método de

rk2, determinando agora não apenas mas também e . Estude parâmetros como o

tempo de queda em função da altura inicial e a a velocidade ao atingir o solo em função da altura de

queda, para diferentes valores de .

5.1.Introdução

Este problema pede-nos para estudar a variação dos tempos e velocidades de queda comparativamente à

altura de queda (v(t), z(t), v(z)), relativamente ao exercício 4. Para podermos criar gráficos comparativos

temos de dar valores a e a . Assim sendo, vamos estudar cada gráfico mantendo um parâmetro

constante e variando o outro.

5.2.Construção e análise dos gráficos:

O primeiro gráfico é o gráfico da posição em função do tempo para diferentes valores de e

para h constante.

27

Desta vez, mantemos constante e variamos a altura inicial h num gráfico da posição em função do

tempo

No terceiro gráfico, variamos , mantendo h constante num gráfico da velocidade em função da

posição.

28

No quarto gráfico, da velocidade em função da posição, mantemos constante e variamos h.

29

5.3.Conclusão

Através da análise dos gráficos posição-tempo, podemos concluir que o tempo de queda é

afectado tanto por como pela altura. afecta o tempo de queda na medida em que, como

visto no problema anterior, quanto maior este for, menor será a velocidade terminal, e quanto

menor for a velocidade mais tempo leva a chegar ao solo. Relativamente à altura, obviamente

que quanto menor esta for, mais perto está o corpo do chão, logo menos tempo levará a

alcançá-lo.

No gráfico velocidade-posição, podemos igualmente observar que quanto menor for maior

será a aceleração do corpo, e consequentemente, maior será a velocidade ao chegar ao solo,

pois como visto anteriormente, maior será a velocidade terminal. No que toca ao gráfico em

que a altura de que o corpo é lançado é variável, podemos verificar que a velocidade com que o

corpo atinge o solo só varia se a altura for tão baixa que o corpo não tenha tempo de atingir a

velocidade terminal, o que acontece quando o corpo é largado de altura de 1m.


Começamos por importar os módulos necessários e definimos as variáveis que serão utilizadas

pelo programa. h0 e gama são variáveis que serão inseridas durante o decorrer do programa.

A função f é a função que devolve os valores da velocidade e da aceleração em função do

instante anterior e será usada por rk2.

from scipy import*

from pylab import*

from numpy import*

v0=0

h0=float(input("Altura inicial?"))

gama=float(input("Gama?"))

t0=0

tf=10.0

npoints=100

h=tf/npoints

x0=array([h0,v0])

yterra=array([[0.,0.],[tf,0.]])

def f(x):

y=zeros(2,float)

y[0]=x[1]

y[1]=gama*x[1]**2-1

return y

30

As funções rk2step e rk2 já são conhecidas dos problemas anteriores.

Por fim, definimos a função plot que nos desenha os três gráficos:

def rk2step(func,x0,h):

k0=func(x0)*h


def rk2(func,t0,x0,h,npoints):

traj=zeros((npoints,len(x0)+1),float)

for i in range(npoints):

traj[i,0]=t0

traj[i,1:len(x0)+1]=x0

x0=rk2step(f,x0,h)

t0+=h

return traj

figure(1)

plot(dados[:,0],dados[:,2])

xlabel('tempo')


title('Queda de um corpo na atmosfera')

figure(2)

plot(dados[:,0],dados[:,1],'r',yterra[:,0],yterra[:,1],'

b')

xlabel('tempo')

ylabel('posicao')


figure(3)

plot(dados[:,1],dados[:,2])

xlabel('posicao')



show()

31

6. Problema 6:

Compare o integrador rk2 com odeint para o problema do oscilador amortecido. Calcule

as diferenças entre as soluções dadas por cada um destes módulos para diferentes valores do

passo de integração usado em rk2.

6.1.Introdução:

Neste exercício é-nos pedidos para comparar os resultados obtidos por dois integradores:

odeint (integrador existente no módulo integrate.scipy) e rk2 (Método de Runge-Kutta de

ordem 2, integrador fornecido nas aulas práticas).

Para compararmos estes dois integradores iremos estudar o problema do oscilador

amortecido que está no texto-guia.

De modo a resolver esta equação numericamente, tal como vimos no texto-guia, reduzimos a

equação a um sistema de duas equações de primeira ordem:

{

Se medirmos o tempo em unidades de √

, ficamos com:

√

√

E obtemos as seguintes expressões:

{

, onde

√

Implementaremos esta função no python para ser integrada por rk2.

def func(x):

y=zeros(2,float)

y[0]=x[1]

y[1]=-x[0]-0.1*x[1]

return y

32

Uma expressão semelhante será usada pelo integrador odeint:

É de esperar que, por usarmos integradores diferentes, as funções terão argumentos diferentes.

No entanto, como se referem às mesmas condições iniciais para a mesma equação, os

resultados serão os mesmos.

6.2.Comparação dos integradores:

A primeira grande diferença entre os dois integradores é que o integrador odeint usa uma

lista de tempos e devolve o valor da função em cada ponto dessa lista, enquanto que o

integrador rk2 utiliza um tempo inicial, e através de um valor de h (o tamanho do passo),

calcula o valor da função no instante , procedendo da mesma forma nos pontos seguintes

até atingir o número de iterações desejado. De modo a poder comparar os dois integradores,

faremos com que odeint aceite uma lista de tempos equivalente aos utilizados por rk2.

Tomaremos os valores de t0=0 e tf=40:

Para o odeint iremos utilizar e para o método de rk2 utilizaremos

h={0.02,0.05,0.1,0.2,0.5} na forma de uma lista, hrk2.

De modo a criar os gráficos, precisamos de armazenar os resultados em tabelas:

Esta função permite-nos criar tabelas de pontos com o h que quisermos da lista de valores de h

que inserirmos. Note-se que o número de pontos varia com o h escolhido, para que os gráficos

tenham o mesmo tamanho.

def f(x,t):

return (x[1],-x[0]-0.1*x[1])

t=arange(t0,tf,h)

hrk2=[0.02,0.05,0.1,0.2,0.5]

def chooseh(i):

matriz=rk2(func,t0,x0,hrk2[i],int(tf/hrk2[i]))

return matriz

dados1=chooseh(0)

dados2=chooseh(1)

...

33

6.3.Comparação e análise dos gráficos:

Agora que já temos uma tabela com os pontos das funções para cada um dos cinco valores

diferentes de h para rk2 e uma tabela com os pontos obtidos por odeint, já podemos

comparar os gráficos, com recurso à função plot:

Neste gráfico já podemos ver consideráveis diferenças quando o tamanho do passo de rk2 é

0.5. De modo a poder comparar melhor os outros valores, fazemos um “zoom” do gráfico:

figure(1)

plot(t,odinty[:,0],'r',dados1[:,0],dados1[:,1],'b',dados

2[:,0],dados2[:,1],'y--',dados3[:,0],dados3[:,1],'g--

',dados4[:,0],dados4[:,1],'r--

',dados5[:,0],dados5[:,1],'b--')

xlabel("tempo")

ylabel("amplitude")

title("Comparacao dos integradores rk2 e odeint")

legend(('odeint h=0.2','rk2 h=0.02','rk2 h=0.05','rk2

h=0.1','rk2 h=0.2','rk2 h=0.5'))

show()

34

Aqui já é possível ver que os vários gráficos começam a divergir. No entanto, rk2 h=0.02 e

odeint h=0.2 permanecem muito próximos. Seria de esperar que odeint h=0.2, tendo

um passo dez vezes maior do que rk2 h=0.02, tivesse um erro maior. Normalmente um

passo maior implica uma pior aproximação. No entanto, isto não se verifica, o que mostra que

odeint é um integrador eficiente.

De modo a melhor comparar os vários gráficos, definimos uma nova função, com o objectivo

de calcular a diferença entre odeint e rk2 para cada valor de h na lista inserida.

Daqui retirámos h=0.5, porque é demasiado divergente e porque a função só aceita tabelas

cuja quantidade de pontos é um múltiplo do número de pontos de odinty.

def diferenca(dados):

dif=zeros((len(odinty),2),float)

c=len(dados)/len(odinty)

for i in range(len(odinty)):

dif[i,0]=dados[c*i,0]

dif[i,1]=dados[c*i,1]-odinty[i,0]

return dif

dif4=diferenca(dados4)




35

Assim, obtemos o seguinte gráfico:

Agora podemos ver claramente a diferença entre os dois integradores para cada valor de h

utilizado em rk2. Quanto mais pequeno for o h, menor será a diferença entre os dois

integradores, de tal modo que a diferença é quase imperceptível para h=0.02. Só se

aproximarmos a imagem é que podemos ver essa diferença:

figure(2)

plot(dif1[:,0],dif1[:,1],'b',dif2[:,0],dif2[:,1],'y',dif

3[:,0],dif3[:,1],'g',dif4[:,0],dif4[:,1],'r')

xlabel("tempo")

ylabel("diferenca")

title("Diferenca entre os integradores rk2 e odeint")

legend(('h=0.02','h=0.05','h=0.1','h=0.5'))

show()

36

6.4.Conclusão

Da análise dos gráficos, podemos concluir que os integradores odeint com h=0.2 e rk2

com h=0.02 são os mais precisos. No entanto, há uma diferença fundamental entre os dois: o

integrador odeint, por ter um passo maior, precisa de menos pontos para desenhar o gráfico,

mas mantém uma grande precisão. rk2, apesar de fazer um gráfico mais “suave”, precisa de

dez vezes mais pontos para obter a mesma precisão que odeint. Isto implica um tempo de

processamento muito maior, o que torna o integrador odeint mais eficiente.


Depois de importar os módulos que serão necessários, definimos as variáveis que farão parte do

nosso código: t0 é o instante inicial, x0 é a posição inicial, h é o tamanho do passo de odeint,

hrk2 é a lista dos valores de h que serão introduzidos em rk2 e tf é o instante final.

pontos garante que os gráficos começam e terminam nos mesmos instantes.

A função f é a função que devolve os valores da velocidade e da aceleração em função do

instante anterior. t dá-nos a lista de tempos a ser utilizada por

odinty=scipy.integrate.odeint(f,x0,t).

Em seguida, definimos a função func que actua como a função f mas com um formato que possa ser

utilizado por rk2.

Definimos a função rk2step e rk2 tal como já foi feito nos outros problemas:

from scipy import*

from pylab import*

import scipy.integrate

t0=0

x0=array([1.,0.])

h=0.2

hrk2=[0.02,0.05,0.1,0.2,0.5]

tf=40

pontos=tf/h

def f(x,t):

return (x[1],-x[0]-0.1*x[1])

t=arange(t0,tf,h)

odinty=scipy.integrate.odeint(f,x0,t)

def func(x):

y=zeros(2,float)

y[0]=x[1]

y[1]=-x[0]-0.1*x[1]

return y

37

A função choose vai-nos devolver uma tabela com os valores das funções integradas por

rk2, tomando como argumento a posição dos elementos na lista hrk2. Depois, definimos as

cinco tabelas correspondentes às funções para cada h.

A função diferenca vai fazer a diferença entre odeint e rk2 para cada ponto. Em

seguida, definimos quatro tabelas com esses valores.

def rk2step(func,x0,h):

k0= func(x0)*h


def rk2(func,t0,x0,h,npoints):

traj=zeros((npoints,len(x0)+1),float)

for i in range(npoints):

traj[i,0]=t0


x0=rk2step(func,x0,h)

t0+=h

return traj

def chooseh(i):

matriz=rk2(func,t0,x0,hrk2[i],int(tf/hrk2[i]))

return matriz

dados1=chooseh(0)

dados2=chooseh(1)

dados3=chooseh(2)

dados4=chooseh(3)

dados5=chooseh(4)

def diferenca(dados):

dif=zeros((len(odinty),2),float)

c=len(dados)/len(odinty)

for i in range(len(odinty)):

dif[i,0]=dados[c*i,0]

dif[i,1]=dados[c*i,1]-odinty[i,0]

return dif





38

Finalmente, usamos a função plot para construir os dois gráficos:

figure(1)

plot(t,odinty[:,0],'r',dados1[:,0],dados1[:,1],'b',dados

2[:,0],dados2[:,1],'y--',dados3[:,0],dados3[:,1],'g--

',dados4[:,0],dados4[:,1],'r--

',dados5[:,0],dados5[:,1],'b--')

xlabel("tempo")

ylabel("amplitude")

title("Comparacao dos integradores rk2 e odeint")

legend(('odeint h=0.2','rk2 h=0.02','rk2 h=0.05','rk2

h=0.1','rk2 h=0.2','rk2 h=0.5'))

figure(2)

plot(dif1[:,0],dif1[:,1],'b',dif2[:,0],dif2[:,1],'y',dif

3[:,0],dif3[:,1],'g',dif4[:,0],dif4[:,1],'r')

xlabel("tempo")

ylabel("diferenca")

title("Diferenca entre os integradores rk2 e odeint")

legend(('h=0.02','h=0.05','h=0.1','h=0.5'))

show()

39

7. Problema 7:

O integrador odeint tem uma maneira mais económica que rk2 de agrupar os seus

argumentos: em vez de to, h e npoints , aceita uma lista de tempos como argumento. Refaça o

código de rk2 para ter a mesma estrutura que odeint.

7.1.Introdução:

Neste problema é-nos pedido para alterar a forma dos argumentos do integrador rk2 de forma a

ficar mais económico e para que fique com a estrutura do integrador odeint. Ou seja, neste

problema vamos pôr o integrador rk2 a funcionar com uma lista de valores t como argumento

ao invés dos argumentos to, h e npoints em que to era o tempo inicial, h o passo de integração e

npoints o número de pontos em que queríamos calcular a solução.

7.2.Como passar os argumentos para a lista de valores t?

Vamos usar uma lista de tempos nos quais queremos calcular a solução da equação diferencial.

No integrador rk2 o primeiro tempo considerado era to aquele em que sabíamos a

solução exacta da equação (a condição inicial que consideramos). Assim, o nosso primeiro

valor da lista t será to.

Para isto teremos de transformar o que se faz facilmente, bastando para isso

recordar ou seja:

{

Assim temos que o primeiro h é sempre to e os seguintes são calculados pelas diferenças dos

dois tempos anteriores, estes valores de h vão ser armazenadas na tabela a que chamaremos, t2.

Basta agora alterar o código para uma indexação correcta, em que o primeiro ponto é o

ponto “solução” e os seguintes são aqueles em que queremos saber a solução, assim basta

indexar estes valores por ordem na tabela.

def rk2(func,x0,t):

t2=[]

for i in range(len(t)):

if i==0:

t2.append(t[i])

else:

t2.append(t[i]-t[i-1])

traj=zeros((len(t),len(x0)+1),float)


h=t2[i]


traj[i,0]=t[i]


return traj

40

7.3.Testando o código:

De modo a testar este código, iremos utilizar a equação diferencial do Problema 1. Desta vez,

em vez de usarmos valores do tempo separados por um h, como no Problema 1, usamos uma

tabela de tempos, o que nos dá mais liberdade de escolha. Definimos ainda a função exct

com a solução exacta para a equação como forma de comparação.

A função plot devolve-nos o seguinte gráfico:

def f(v):

return -v

t=[0.0,0.1,0.2,0.21,0.22,0.221,0.222,0.223,0.225,0.3,0.4

,0.6,0.9,0.95,0.97,0.98,1.]

x0=array([[1.0]])

exct=zeros((len(t),2),float)


exct[i,0]=t[i]

exct[i,1]=exp(-t[i])

rk2alt=rk2(f,x0,t)

figure(1)

plot(exct[:,0],exct[:,1],'b',rk2alt[:,0],rk2alt[:,1],'r-

o')

xlabel('tempo')

ylabel('altura')

title('rk2alterado')

legend(('Exacto','rk2alterado'))

show()

41

7.4. Conclusão:

Conseguimos alterar o código de rk2 com sucesso e agora esta função trabalha com os mesmos

argumentos da função odeint.

Podemos agora dar o ponto solução em qualquer ponto desde que saibamos o primeiro t (o t do

ponto “solução”). Note-se contudo, que este método depende do valor de h para fazer uma boa

aproximação. Por isso, não convém escolher pontos muito separados uns dos outros.

7.5.Análise do código:

Começamos por importar os módulos scipy, pylab e numpy, e definimos rk2step, o

passo do Método de Runge-Kutta de ordem 2. Em seguida, definimos rk2 que começa por

criar uma nova lista com as diferenças entre os valores da lista de tempos inserida. Em seguida,

rk2 vai usar essa nova lista para criar uma tabela com os valores da função nesses pontos.

from scipy import *

from pylab import *

from numpy import *

def rk2step (func,x0,h):

k0=func(x0)*h


def rk2(func,x0,t):

t2=[]


if i==0:

t2.append(t[i])

else:

t2.append(t[i]-t[i-1])

traj=zeros((len(t),len(x0)+1),float)


h=t2[i]


traj[i,0]=t[i]


return traj

42

Em seguida, para testar, inserimos uma lista t com os tempos onde queremos calcular a

solução da equação o valor inicial x0 da função, e definimos a função f que é a derivada da

função.

Definimos ainda a solução exacta para a equação:

Finalmente, utilizamos a função plot para desenhar o gráfico:

t=[0.0,0.1,0.2,0.21,0.22,0.221,0.222,0.223,0.225,0.3,0

.4,0.6,0.9,0.95,0.97,0.98,1.]

x0=array([[1.0]])

def f(v):

return -v

rk2alt=rk2(f,x0,t)

exct=zeros((len(t),2),float)


exct[i,0]=t[i]

exct[i,1]=exp(-t[i])

figure(1)

plot(exct[:,0],exct[:,1],'b',rk2alt[:,0],rk2alt[:,1],'

r-o')

xlabel('tempo')

ylabel('altura')

title('rk2alterado')

legend(('Exacto','rk2alterado'))

show ( )

Documents

Trabalho de Práticas Computacional de Mecânica