Assistance with adaptive time step for Dormand Prince Method

Click For Summary
SUMMARY

The discussion focuses on implementing an adaptive time step for the Dormand-Prince method in a MATLAB program designed to simulate solar system dynamics. The user successfully created a Runge-Kutta 4th order integrator but encounters issues with excessively small time steps, specifically a dt value of 4.2431e-06, leading to prolonged computation times. The user is advised to ensure the correct calculation of position and velocity differences (difp and difv) and to utilize the minimum value of these differences to optimize the time step effectively.

PREREQUISITES
  • Understanding of the Dormand-Prince method for numerical integration
  • Familiarity with MATLAB programming and syntax
  • Knowledge of Newton's Law for calculating acceleration
  • Experience with adaptive time-stepping techniques in numerical methods
NEXT STEPS
  • Explore MATLAB's ode45 solver for insights on adaptive time-stepping
  • Research error estimation techniques in numerical integration
  • Learn about optimizing performance in MATLAB for large-scale simulations
  • Investigate common pitfalls in implementing the Dormand-Prince method
USEFUL FOR

Students and researchers in computational physics, MATLAB programmers, and anyone interested in numerical methods for simulating dynamic systems.

jbrisby
Messages
5
Reaction score
0

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))
 
Physics news on Phys.org
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!
 
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.