TRABALHO 6 MÉTODO DE RK4 PARA SISTEMAS COM RK3

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

clc
clear
//COMPARANDO RK4 COM RK3
//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 para ambos
solrk4(1,1)=10
solrk4(2:30,1)=0
solrk3(1,1)=10
solrk3(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 rk3
while t(i)<=tf
 //calculando com rk3
 kr1=f(solrk3(1:30,i))
 kr2=f(solrk3(1:30,i)+dt*kr1/2)
 kr3=f(solrk3(1:30,i)-dt*kr1+2*dt*kr2)
 solrk3(1:30,i+1)=solrk3(1:30,i)+dt/6*(kr1+4*kr2+kr3)
 //calculando com rk4
 k1=f(solrk4(1:30,i))
 k2=f(solrk4(1:30,i)+k1*dt/2)
 k3=f(solrk4(1:30,i)+k2*dt/2)
 k4=f(solrk4(1:30,i)+k3*dt)
 solrk4(1:30,i+1)=solrk4(1:30,i)+dt/6*(k1+2*k2+2*k3+k4)
 //erro comparado entre os dois metodos
 erro=max(abs(solrk4(1:30,i+1)-solrk3(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(solrk4)