Runge Kutta fehlberg not going through full simulation
3 views (last 30 days)
Show older comments
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 Comments
Accepted Answer
Fei Deng
on 23 Sep 2016
Edited: Fei Deng
on 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:
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!