Runge Kutta fehlberg not going through full simulation
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi guys, Hope you can help me with an issue I'm currently having with a RKF45 simulation.
My code is copied below. Basically, I've got a 4th order Runge Kutta which works fine and gives me 86400 predictions to an ODE. I want the Runge Kutta Fehlberg to do the same (hopefully more accurately though) but it only gives me 2705 predictions. I would like to have 86400. Any ideas why? I can't see what's missing... Cheers
if true
f= @(R_moon) (-G*(moon_mass+earth_mass)*R_moon)/norm(R_moon)^3;
R1=moon_position_vector(2,:);
V1=data_moon(2,8:10);
t=1; tfin=24*3600;
h=1; hmin=h/64; hmax=64*h;
j=1;
epsilon=0.0000001;
max1=200;
while (t<=tfin)
k1 = h*f(R1);
k2 = h*f(R1+(1/4)*k1);
k3 = h*f(R1+(3/32)*k1+(9/32)*k2);
k4 = h*f(R1+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3);
k5 = h*f(R1+(438/216)*k1-8*k2+(3680/513)*k3-(845/4104)*k4);
k6 = h*f(R1-(8/27)*k1+2*k2-(3544/2565)*k3+(1859/4104)*k4-(11/40)*k5);
err=max(abs(((1/360)*k1-(128/4275)*k3-(2197/75240)*k4+(1/50)*k5+(2/55)*k6)));
Vnew= V1+(25/216)*k1+(1408/2565)*k3+(2197/4101)*k4-(1/5)*k5;
Vnew2 = V1 + 16*k1/135 + 6656*k3/12825 + 28561*k4/56430 - 9*k5/50 + 2*k6/55;
R=max(abs(Vnew-Vnew2));
delta=0.84*(epsilon*h/R)^(0.25);
if ((err<epsilon)|(h<2*hmin))
V(j,:)=Vnew2;
R1=R1+Vnew*h;
Rnew(j,:)=R1;
end
t=t+h;
time(j) = t;
j=j+1;
if ((delta<0.75)&(h>2*hmin)), h = h/2; end
if ((delta>1.50)&(2*h<hmax)), h = 2*h; end
end
end
0 commentaires
Réponse acceptée
Fei Deng
le 23 Sep 2016
Modifié(e) : Fei Deng
le 23 Sep 2016
Hi Joe,
RKF45 method allows for an adaptive step size to be determined automatically, which improves efficiency/reduce steps. I suggest you first take a look at the way this method works.
Moreover, there is a function to realize RKF45 method in MATLAB File Exchange:
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!