Runge Kutta fehlberg not going through full simulation

Question

Joe O'Leary 2016-9-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/303404-runge-kutta-fehlberg-not-going-through-full-simulation

评论： Joe O'Leary 2016-9-24

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Fei Deng 2016-9-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/303404-runge-kutta-fehlberg-not-going-through-full-simulation#answer_235855

编辑：Fei Deng 2016-9-23

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:

https://www.mathworks.com/matlabcentral/fileexchange/2181-numerical-methods-using-matlab-2e/content/edition2/matlab/chap_9/rkf45.m

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Joe O'Leary 2016-9-24

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.

请先登录，再进行评论。