Assistance with adaptive time step for Dormand Prince Method

[jbrisby]

Homework Statement

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.

Homework Equations

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

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))

 

[AlephZero]

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!

 

[jbrisby]

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.