Daniel Carrera <email@example.com> writes: >Hello, > >I'm making a personal library of ODE integrators. I want to have three >integrators: > >1) A traditional 4th-order Runge-Kutta integrator. Simplest method, >easiest to get right. Ideally, I'd like to add an adaptive step size. > >2) A 4th+5th order adaptive step size integrator of the Runge-Kutta >family. Like the Dormand-Prince or similar. > >3) A Bulirsch-Stoer method with Richardson extrapolation. > > >I have two questions pertaining the first two integrators: > >1) I already have the basic 4th-order RK method implemented. I think >that the simplest reasonable way to give it an adaptive step size is >to use the mid-point method, like this: > >! >! Variables needed for 4th-order Runge-Kutta >! >k1 = dY(t, Y) >k2 = dY(t + h/2, Y + k1*h/2) >k3 = dY(t + h/2, Y + k2*h/2) >k4 = dY(t + h , Y + k3*h) > >! >! The second also gives me the mid-point method. >! >Y_rk4 = Y + (k1 + 2*k2 + 2*k3 + k4) * h/6 >Y_mid = Y + k2 * h > >! >! Error estimate and new step size. >! >error = abs(Y_rk4 - Y_mid) / ( abs(Y_rk4) * rtol ) > >h = h * (rtol / error) ** (1/2) > >if (error > rtol) then > ... go back and try again ... >end if > > >Does this look reasonable?
hmmm, you embed a second order scheme in the fourth order scheme, hence you estimate the error in the second order scheme and integrate using its order (2) and stbility bound (which is smaller than the one of RK4.) This will be inefficient . >
some safeguards should be added: reduction not too rapid, a smallest step size introduced , you also will need to increase the stepsize (with safeguards) , a max stepsize must be introduced,
> >2) My second question is shorter. There are several 4th+5th order >methods of the Runge-Kutta family: > >* Dormand-Prince >* Fehlberg >* Cash-Karp > >Which one should I pick? > >Thanks for the help. > >Daniel. Dormand-Prince did minimize the error constant in the 5th order scheme, Fehlberg is inferior, I don't know the stability bound of Cash-Karp, if this is larger than the one of Dormand-Prince, it may be advantageous. Hence: it depends on the stiffness you may expect from your systems. (There is also a good 4/5th order schem by Shampine, quite good ..) hth peter