Exercícios de Fisica computacional

Prévia do material em texto

Universidade Federal do Oeste da Bahia – UFOB 
Centro de Ciências Exatas e das Tecnologias – CCET 
Exercícios de Física Computacional – Semestre Letivo 2016.1 
 
FÍSICA COMPUTACIONAL 
Prof.: Antonio César do Prado Rosa Junior 
Aluno: Márcio Moreira Lima 
Curso: Física 
 
Com base nos algoritmo já feito nos exemplos do livro Métodos Computacional da Física, Claudio Sherer. 2a 
edição. Utilizamos o Scilab para desenvolver algoritmos para resolver os exercícios proposto no livro. 
 
 
Exemplo 2.1 
 
Implementar a solução numérica de. 
 
dx/dt = x x(0)=1 
Com ambas as formas e comparar com a solução exata. 
 Solução: A solução é feita de duas formas distintas de 
derivada numérica, à direita e a centrada. 
Pela direita temos que: 
xj+1 = xj + dt*xj 
xj+1 = xj (1+ dt) 
Onde xj+1 é o valor futuro e xj é valor atual. 
Para a derivada numérica centrada temos. 
xj+1 = xj-1 + dt*xj 
Temos que fazer algumas correções na equação 
devidas o valor de x1 que não temos para calcular x2. 
Então podemos reescrever como sendo. 
xj+1 = xj-1 + dt*(x0 (1+dt)) 
A solução exata é: 
x = exp(t) 
 
 
 
A seguir esta o algoritmo implementado em Scilab. 
 
derivadas.sce 
 
tmax=5 
dt=0.01 
 
nt= floor (tmax/dt)+1 
t=linspace (0,tmax,nt); 
x=zeros(1,nt); 
x(1)=1; 
umdt=1+dt; 
h=2*dt; 
for n=2:nt; 
x(n)= umdt*x(n-1); 
end 
y(1)=1; 
y(2)=umdt*x(1); 
for n=3:nt; 
y(n)=y(n-2)+h*y(n-1); 
end 
z=exp(t); 
xset("window", 1); 
clf; 
plot(t,x,t,y,t,z,"--") 
legend ("a direita","centrada", "exata",2) ; 
 
 Usando a condições iniciais tmax=5 e a variação do 
passo dt=1 entramos o seguinte gráfico. 
 Figura 1- Gráfico de (t , x) obtido a partir dt =0.1 e tmax =1. 
 Figura 2: Gráfico de (t , x) obtido a partir de dt = 0.01 e tmax =5 
 
Podemos observar na figura 1 e 2 que quando aumentamos o 
dt à derivada numérica a direito converge para a solução 
exata e podemos perceber que a derivada centrada é mais 
precisa nesta situação. 
 
Exercício 2.1: Modifique o programa para resolver a 
equação. 
dx/dt = a * x² 
Com x(0) = 1, e execute-o com diferentes valores da 
constante a, positivos e negativos, comparando com a 
relação exata x = 1 /(1- a *t). 
 
Solução: com base no exemplo 2.1 vamos construir nosso 
algoritmo. 
 
 
 
 
Universidade Federal do Oeste da Bahia – UFOB 
Centro de Ciências Exatas e das Tecnologias – CCET 
Exercícios de Física Computacional – Semestre Letivo 2016.1 
 
tmax = ?. 
dt = ? s 
a= ? 
nt= floor (tmax/dt)+1 
t=linspace (0,tmax,nt); 
x=zeros(1,nt); 
x(1)=1; 
h=2*a*dt; 
for n=2:nt; 
 x(n)=x(n-1)*(1+a*dt*x(n-1)); 
end 
y(1)=1; 
y(2)= x(1)*(1+a*dt*x(1)); 
for n=3:nt; 
 y(n)=y(n-2)+h*y(n-1)^2; 
end 
z=(1-a*t).\1; 
xset("window", 1); 
clf; 
plot(t,x,t,y,t,z,"--") 
legend ("a direita","centrada", "exata",2) 
xtitle ("Gráfico para a = 0.1 e tmax=5") 
 
Para a >1 a solução exata diverge. Vamos iniciar 
procedimento numérico para t = 0 e a*t < 1. Vamos ver os 
gráficos abaixo para diferentes valores de tmax e a. 
 Figura 3- Gráfico de (t , x) para os valores de a = 0.1 e tmax = 5
 Figura 4- Gráfico obtido com valores de a = 0.1 e tmax = 9.9. 
 
Como podemos observar na figura 3 e 4 só obtemos solução 
para a > 0 se a*t < 1. Para tmax = 10 a equação diverge e 
para tmax = 9.99 só conseguimos obter solução para valores 
de a <= 0.1. 
 
Para a < 0 obtemos os seguintes gráficos: 
 
 
 Figura5- Gráfico de (t , x) a partir de a = - 0.1 e tmax = 5 
 
 Figura6 - Gráfico de (t , x) para valores de a = - 0.1 e tmax = 10. 
 
 
 Figura7- Gráfico para valores de a = -1 e tmax = 10. 
 
Para valores de a < 0 a derivada centrada apresenta 
comportamento estranho se o |𝑎 ∗ 𝑡| não for muito pequeno 
Universidade Federal do Oeste da Bahia – UFOB 
Centro de Ciências Exatas e das Tecnologias – CCET 
Exercícios de Física Computacional – Semestre Letivo 2016.1 
 
para a = -1 e tmax = 10 a derivada centrada é diferente do 
esperado pela equação diferencial que se pretende encontrar. 
Podemos perceber nesse caso que a derivada centrada não é 
indicada para resolver o problema. 
 
Exemplo 2.2: seja x = (x,y) f = (-ay,bx). A equação 
equivalente a equação xj+1 = xj + dt*xj é o sistema. 
xj+1 = xj – a*dt*yj 
 yj+1 = yj + b*dt*xj 
Implementando o exemplo acima temos e podemos 
definir x como vetor coluna, os passos de integração com 
auxílio da matriz. 
 
Sistema.sce. 
a=1 
b=1 
x0=1 
y0=0 
tmax= 12 
dt=0.01 
 
X0=[x0;y0] 
nt=floor(tmax/dt)+1; 
t=linspace(0,tmax,nt); 
X=zeros(2,nt); 
X(:,1)= X0; 
A=[0,-a;b,0] 
for j=1:nt-1; 
 X(:,j+1)=X(:,j)+dt*A*X(:,j); 
end 
x=X(1,:); 
y=X(2,:); 
clf; 
plot(t,x,t,y); 
legend("x(t)","y(t)"); 
 
o gráfico obtido é: 
 Figura 8- Gráfico de (t,x; t,v) gerado a partir função vetorial. 
 
 
 
 
 
 
Exercícios 2.2 
Modifique o programa sistema.sc, usando para f uma função 
não linear x e y. por exemplo, fx = -ay³, fy = bx³. 
 Solução: 
xj+1 = xj – a*dt*yj³ 
 yj+1 = yj + b*dt*xj³ 
o algoritmo é 
a=1 
b=1 
X0=1 
Y0=0 
tmax= 12 
dt=0.01 
 
nt=floor(tmax/dt)+1; 
t=linspace(0,tmax,nt); 
X(1)= X0; 
Y(1)=Y0; 
for j=1:nt-1; 
 X(j+1)=X(j)-dt*a*Y(j)^3; 
 Y(j+1)=Y(j)+dt*b*X(j)^3; 
end 
plot(t,X,t,Y); 
legend("x(t)","y(t)"); 
 
O gráfico obtido é: 
 Figura 9- Gráfico de (t,x; t,v) da função vetorial. 
 
Exemplo 3.3 
Considere um pendulo de comprimento l=1. A equação de 
movimento para o ângulo Ө, entre a barra pendulo e a 
vertical é: 
Ӫ = -g sin(Ө) 
Solução: a seguir programaremos o algoritmo no Scilab para 
resolver numericamente o exemplo: 
 
clear 
close 
clc 
teta0= %pi/4 
v0=5 
tmax= 30 
dt=0.01 
ddt=2*dt 
dt2=dt^2 
nt= floor(tmax/dt)+1; 
t=linspace(0,tmax,nt); 
teta(1)=teta0 
v(1)= v0 
teta(2)= v(1)*dt+teta(1) ; 
v(2)= (teta(2)-teta(1))/dt 
for j=2:nt-1; 
 f(j)= -10*sin(teta(j)); 
Universidade Federal do Oeste da Bahia – UFOB 
Centro de Ciências Exatas e das Tecnologias – CCET 
Exercícios de Física Computacional – Semestre Letivo 2016.1 
 
 teta(j+1)=2*teta(j)-teta(j-1)+((f(j)*dt2)); 
 v(j+1)=(teta(j+1)-teta(j-1))/ddt; 
end 
 
xset("window",1); 
clf; 
plot(t,teta,t,v); 
legend("Verlet","velocidade",2); 
 
Podemos ver abaixo o gráfico do pendulo simples obtido a 
partir do algoritmo. 
 
 Figura 10- Gráfico de (t , x; t , v) obtido para Ө= 𝟏. 𝟓 𝒆 𝒕𝒎𝒂𝒙 = 𝟏𝟎 𝒆 𝒗𝟎 =
𝟎. 
 
 
Exemplo 4. 
Generalizar o algoritmo do pendulo simples para o 
caso do oscilador massa-mola de m=1. Pela equação 
dada: 
𝑑²𝑥
𝑑𝑡²
= −𝑘𝑥 
O algoritmo em Scilab é: 
 
clear 
close 
clc 
x0=0 
v0=5 
k= 0.91 
tmax= 30 
dt=0.01 
ddt=2*dt 
dt2=dt^2 
nt= floor(tmax/dt)+1; 
t=linspace(0,tmax,nt); 
x(1)=x0 
v(1)= v0 
x(2)= v(1)*dt+x(1) ; 
v(2)= (x(2)-x(1))/dt 
for j=2:nt-1; 
 f(j)= -k*x(j) 
 x(j+1)=2*x(j)-x(j-1)+((f(j)*dt2)); 
 v(j+1)=(x(j+1)-x(j-1))/ddt; 
end 
 
xset("window",1); 
clf; 
plot(t,x,t,v); 
legend("trajetória Verlet ","velocidade Valet",-1); 
 
Segue abaixo o gráfico a partir do algoritmo: 
 
 Figura 11- Gráfico de (t,x; t,v) para caso MHS. 
 
Exemplo 5. 
Vamos tratar do caso de uma função dependente da 
velocidade que o casa de uma partícula sobre ação de 
uma forçade atrito viscoso descrito pela seguinte 
equação: 
𝑑²𝑥
𝑑𝑡²
= −𝛾𝑣 
Solução: 
O usando método de Verlet pode fazer o algoritmo em 
Scilab como segue abaixo: 
 
clear 
close 
clc 
x0=2 
v0=10 
k= 0.25 
tmax= 6 
dt=0.01 
ddt=2*dt 
dt2=dt^2 
nt= floor(tmax/dt)+1; 
t=linspace(0,tmax,nt); 
x(1)=x0 
v(1)= v0; 
x(2)= v(1)*dt+x(1) ; 
v(2)= (x(2)-x(1))/dt 
for j=2:nt-1; 
 f(j)= -k*(v(j)); 
 x(j+1)=2*x(j)-x(j-1)+((f(j)*dt2)); 
 v(j+1)=(x(j+1)-x(j))/dt; 
end 
z=v0*exp(-k*t); 
xset("window",1); 
clf; 
plot(t,x,t,v,t,z,"--"); 
legend("Verlet","velocidade","Quase exata",-1); 
Universidade Federal do Oeste da Bahia – UFOB 
Centro de Ciências Exatas e das Tecnologias – CCET 
Exercícios de Física Computacional – Semestre Letivo 2016.1 
 
 
Podemos observar os gráficos a baixo gerado pelo 
programa. 
 
 Figura 12- Gráfico de (t , x ; t , v) para o atrito viscoso. 
 
 
 Figura 13- Gráfico de (t ,x ; t , v) para o atrito viscoso com 𝜸 = 𝟏. 
 
Exemplo 6. 
No caso mais simples onde temos a seguinte equação: 
𝑑²𝑥
𝑑𝑡²
= 𝑘 
 
Solução: O usando método de Verlet pode fazer o 
algoritmo em Scilab como segue abaixo: 
clear 
close 
clc 
x0=0 
v0=12 
k= 0.45 
tmax= 5 
dt=0.1 
ddt=2*dt 
dt2=dt^2 
nt= floor(tmax/dt)+1; 
t=linspace(0,tmax,nt); 
x(1)=x0 
v(1)= v0 
x(2)= v(1)*dt+x(1) ; 
v(2)= (x(2)-x(1))/dt 
for j=2:nt-1; 
 f(j)= k 
 x(j+1)=2*x(j)-x(j-1)+((f(j)*dt2)); 
 v(j+1)=(x(j+1)-x(j))/dt; 
end 
z=v0+k*t; 
xset("window",1); 
clf; 
plot( t,x,t,v,t,z,"--"); 
legend("Posição Verlet","velocidade Verlet","Quase exata",2); 
 
A seguir podemos observar o gráfico: 
 
 Figura 14- Gráfico de (t , x; t , v) para uma constate.