Azzera filtri
Azzera filtri

Runge Kutta fehlberg not going through full simulation

7 visualizzazioni (ultimi 30 giorni)
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

Risposta accettata

Fei Deng
Fei Deng il 23 Set 2016
Modificato: Fei Deng il 23 Set 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:
  1 Commento
Joe O'Leary
Joe O'Leary il 24 Set 2016
Thank you Fei. I realised my naive question quite soon after posting it. I have adapted the RFK45 function for a two body system which is working quite well now. Thanks again.

Accedi per commentare.

Più risposte (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by