Runge Kutta fehlberg not going through full simulation
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Joe O'Leary
am 17 Sep. 2016
Kommentiert: Joe O'Leary
am 24 Sep. 2016
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 Kommentare
Akzeptierte Antwort
Fei Deng
am 23 Sep. 2016
Bearbeitet: Fei Deng
am 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:
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!