TRABALHO 6 MÉTODO DE RK4 PARA SISTEMAS COM EULER

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

clc
clear
//COMPARANDO RK4 COM EULER
//definir função
function fun=f(y) //y é um vetor de 30 componentes y(i)=xi
 fun(1)=0 //fun é um vetor de 30 comp
 for i=2:29
 fun(i)=4*(y(i+1)-2*y(i)+y(i-1))-0.95*(y(i+1)-y(i-1))-0.5*y(i)
 end
 fun(30)=fun(29)
endfunction
//definir condições iniciais
solrk(1,1)=10
solrk(2:30,1)=0
sole(1,1)=10
sole(2:30,1)=0
t(1)=0
//definir parametros para while, dt variável
dt=0.1
tf=1
prec=0.05
i=1
//resolução por rk4 comparado com euler
while t(i)<=tf
 //calculando euler
 ke=f(sole(1:30,i))
 sole(1:30,i+1)=sole(1:30,i)+ke*dt
 //calculando rk4
 k1=f(solrk(1:30,i))
 k2=f(solrk(1:30,i)+k1*dt/2)
 k3=f(solrk(1:30,i)+k2*dt/2)
 k4=f(solrk(1:30,i)+k3*dt)
 solrk(1:30,i+1)=solrk(1:30,i)+dt/6*(k1+2*k2+2*k3+k4)
 //erro comparando os dois metodos
 erro=max(abs(solrk(1:30,i+1)-sole(1:30,i+1)))
 if erro>prec
 //dt=((prec/erro)**0.25)*dt - nao roda
 dt=0.7*dt
 else
 //dt=((prec/erro)**0.2)*dt
 dt=1.3*dt
 t(i+1)=t(i)+dt
 i=i+1
 end
end
disp(solrk)