lista 3 - 2015.2

Prévia do material em texto

Pontificia Universidade Católica do Rio de Janeiro
 Departamento de Engenharia Mecânica
 
 
 
 ENG 1714 
 Métodos Numéricos para Engenharia Mecânica
Aluno: Leonardo Pedreira Pereira
Matricula: 1220832
Turma: 3VB
Professor: Márcio Carvalho
Lista 3
Ex 1) Sistema massa-mola dado por:
Reduzindo a um sistema de EDO com x =x1:
SCI LAB (Metodo Runge Kutta 4 ordem)
function y=F(x2)
 y=x2
endfunction
function y=G(x1, x2, t, w, Fo, m, a, k)
 y=Fo*sin(w*t)/m-k*x1/m-a*abs(x2)*x2/m
endfunction
m=2
Fo=2
a=5
k=5
tf=15
i=1
x1(1)=1
x2(1)=0
t(1)=0
h=input('Passo de tempo: ')
w=input('Frequencia: ')
while t(i)<=tf
 t(i+1)=t(i)+h
 
 f1=F(x2(i))
 g1=G(x1(i),x2(i),t(i),w,Fo,m,a,k)
 f2=F(x2(i)+g1*h/2)
 g2=G(x1(i)+f1*h/2,x2(i)+g1*h/2,t(i)+h/2,w,Fo,m,a,k)
 f3=F(x2(i)+g2*h/2)
 g3=G(x1(i)+f2*h/2,x2(i)+g2*h/2,t(i)+h/2,w,Fo,m,a,k)
 f4=F(x2(i)+g3*h)
 g4=G(x1(i)+f3*h,x2(i)+g3*h,t(i)+h,w,Fo,m,a,k)
 
 x1(i+1)=x1(i)+(h/6)*(f1+2*f2+2*f3+f4)
 x2(i+1)=x2(i)+(h/6)*(g1+2*g2+2*g3+g4)
 
 i=i+1
end
figure
plot(t,x1,t,x2);
legend('deslocamento[m]','velocidade[m/s]')
a)Resultado na tela:
Passo de tempo: 0.01
Frequencia: 0.5
B)Resultado na tela:
Passo de tempo: 0.01
Frequencia: 2
As soluções independem do tamanho passo do tempo.Se colocarmos um passo de tempo maior para as soluçoes, o gráfico será o mesmo porém ficará distorcido. Se usarmos um intervalo de tempo menor, o gráfico ficará encolhido. Um exemplo disso é no gráfico da letra A e B usando h=0.9 
Ex 2) 
Metodo de Euler Implícito: 
Como F nao é uma função linear, resolveremos pelo método de newton.
SCI LAB (Método de Euler Implícito de 1 ordem + Método de Newton)
 function z=F(y, k)
 z=-k*(y)^(1/2)
endfunction
function z=dF(y, k)
 z=-k/(2*(y)^(1/2))
endfunction
function z=newton(b, Dt, k)
 y0=b
 erro=abs(y0-b-Dt*F(y0,k))
 while(erro>1e-4)
 Dy=-(y0-b-Dt*F(y0,k))/(1-Dt*dF(y0,k))
 y0=y0+Dy
 erro=abs(y0-b-Dt*F(y0,k))
 end
 z=y0
endfunction
ci=input('Nivel de agua: ')
Dt=input('Digite o passo de tempo: ')
k=input('Digite o valor de k: ')
t(1)=0
y(1)=ci
yf=0.0001
i=1
while y(i)>=yf
 t(i+1)=t(i)+Dt
 y(i+1)=newton(y(i),Dt,k)
 i=i+1
end
disp('Tempo para esvaziar o tanque: ')
disp(t(i))
Nivel de agua: 3
Digite o passo de tempo: 0.1
Digite o valor de k: 0.1
 Tempo para esvaziar o tanque: 
 34.7