1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Assistance with adaptive time step for Dormand Prince Method

  1. Jul 28, 2014 #1
    1. The problem statement, all variables and given/known data
    I am trying to create a matlab program that propagates the solar system when given initial velocity and position. I have successfully implemented a Runge-Kutta 4th order integrator but would like to see the difference when moving to an adaptive time step method. I think my problem arises when it comes to choosing the optimum time step for the next iteration. When checking the error of the position and using that to calculate the next time step I get a dt that approaches 4.2431e-06 and causes the integration to take exceedingly long periods of time to finish. The difp and difv are matrices of the error for position and velocity. I assume that I'm going to take the minimum value of the position difference because it will require the smaller time step. Any help would be greatly appreciated. Also, if it would help I can elaborate or provide more information if needed.


    2. Relevant equations

    A reference for the Dormand Prince method if needed: http://depa.fquim.unam.mx/amyd/archivero/DormandPrince_19856.pdf

    3. The attempt at a solution
    In the following code my accel() function calculates acceleration using Newton's Law;
    pos and vel are 11x3 matrices for the bodies of the solar system being calculated;
    dt is the initial time step;
    c is a matrix of coefficients for the integration method;
    c = [0 0 0 0 0 0
    1/5 0 0 0 0 0
    3/40 9/40 0 0 0 0
    44/45 -56/15 32/9 0 0 0
    19372/6561 -25360/2187 64448/6561 -212/729 0 0
    9017/3168 -355/33 -46732/5247 49/176 -5103/18656 0
    35/384 500/1113 125/192 -2187/6784 11/84 0
    5179/57600 7571/16695 393/640 -92097/339200 187/2100 1/40];

    k1v = accel(pos)*dt;
    k1p = vel*dt;
    k2v = accel(pos+(k1p.*c(2,1))).*dt;
    k2p = (vel+(k1v.*c(2,1))).*dt;
    k3v = accel(pos+(k1v.*c(3,1)+k2v.*c(3,2))).*dt;
    k3p = (vel+(k1p.*c(3,1)+k2p.*c(3,2))).*dt;
    k4v = accel(pos+(k1v.*c(4,1)+k2v.*c(4,2)+k3v.*c(4,3))).*dt;
    k4p = (vel+(k1p.*c(4,1)+k2p.*c(4,2)+k3p.*c(4,3))).*dt;
    k5v = accel(pos+(k1v.*c(5,1)+k2v.*c(5,2)+k3v.*c(5,3)+k4v.*c(5,4))).*dt;
    k5p = (vel+(k1p.*c(5,1)+k2p.*c(5,2)+k3p.*c(5,3)+k4p.*c(5,4))).*dt;
    k6v = accel(pos+(k1v.*c(6,1)+k2v.*c(6,2)+k3v.*c(6,3)+k4v.*c(6,4)+k5v.*c(6,5))).*dt;
    k6p = (vel+(k1p.*c(6,1)+k2p.*c(6,2)+k3p.*c(6,3)+k4p.*c(6,4)+k5p.*c(6,5))).*dt;
    k7v = accel(pos+(k1v.*c(7,1)+k3v.*c(7,2)+k4v.*c(7,3)+k5v.*c(7,4)+k6v.*c(7,5))).*dt;
    k7p = (vel+(k1p.*c(7,1)+k3p.*c(7,2)+k4p.*c(7,3)+k5p.*c(7,4)+k6p.*c(7,5))).*dt;
    k8v = accel(pos+(k1v.*c(8,1)+k3v.*c(8,2)+k4v.*c(8,3)+k5v.*c(8,4)+k6v.*c(8,5)+k7v.*c(8,6))).*dt;
    k8p = (vel+(k1p.*c(8,1)+k3p.*c(8,2)+k4p.*c(8,3)+k5p.*c(8,4)+k6p.*c(8,5)+k7p.*c(8,6))).*dt;
    difv = abs(k8v-k7v);
    difp = abs(k8p-k7p);
    s=(tol*dt./(2.*difp)).^.2;
    dt = min(min(s.*dt))
     
  2. jcsd
  3. Jul 28, 2014 #2

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    According to http://en.wikipedia.org/wiki/Dormand–Prince_method (referenced in the PDF) this is the same method that is used in the ode45 solver in Matlab.

    If that is correct (and you might want to confirm it in the Matlab documenation) and your code gives different results, something must be wrong somewhere.

    There are lots of opportunities to make a typo in the style of code you attached, but I'm not going to check every character to try and find it!
     
  4. Jul 28, 2014 #3
    Yeah it is the one that MATLAB uses but my prof wants me to write my own as a "learning opportunity"... heh. And no worries I know it's a pain debugging other people's code. I think my problem is something with choosing the correct difference at the end. difv and difp both give an 11x3 matrix of difference values. So I took the minimum value out of all those and tried to use that to make the time step. Thanks for the reply though.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Assistance with adaptive time step for Dormand Prince Method
  1. Assistance Required (Replies: 2)

  2. Please assist! (Replies: 1)

  3. Paid Assistance? (Replies: 2)

Loading...