Skip to content

Commit

Permalink
added higher-order Adams formula.
Browse files Browse the repository at this point in the history
added up to 8th order Adams to oscillation test.
  • Loading branch information
jacobwilliams committed Dec 30, 2015
1 parent 4bdee14 commit 360c25a
Show file tree
Hide file tree
Showing 4 changed files with 310 additions and 38 deletions.
168 changes: 152 additions & 16 deletions src/lib/foodie_integrator_adams_bashforth.f90
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
!< FOODIE integrator: provide an explicit class of Adams-Bashforth multi-step schemes, from 1st to 4rd order accurate.
!< FOODIE integrator: provide an explicit class of Adams-Bashforth multi-step schemes, from 1st to 16th order accurate.
module foodie_integrator_adams_bashforth
!-----------------------------------------------------------------------------------------------------------------------------------
!< FOODIE integrator: provide an explicit class of Adams-Bashforth multi-step schemes, from 1st to 4rd order accurate.
!< FOODIE integrator: provide an explicit class of Adams-Bashforth multi-step schemes, from 1st to 16th order accurate.
!<
!< Considering the following ODE system:
!<
Expand Down Expand Up @@ -54,6 +54,7 @@ module foodie_integrator_adams_bashforth
!< +\frac{37}{24} R(t^{n+1}, U^{n+1}) - \frac{9}{24} R(t^{n}, U^{n}) \right] $$
!<
!<#### Bibliography
!< * J. L. Maury Jr., G. P. Segal, "[Cowell Type Numerical Integration As Applied to Satellite Orbit Computation](http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19690017325.pdf)", NASA-TM-X-63542, X-553-69-46, April 1969.
!<
!-----------------------------------------------------------------------------------------------------------------------------------

Expand All @@ -70,7 +71,7 @@ module foodie_integrator_adams_bashforth

!-----------------------------------------------------------------------------------------------------------------------------------
type :: adams_bashforth_integrator
!< FOODIE integrator: provide an explicit class of Adams-Bashforth multi-step schemes, from 1st to 4rd order accurate.
!< FOODIE integrator: provide an explicit class of Adams-Bashforth multi-step schemes, from 1st to 16th order accurate.
!<
!< @note The integrator must be created or initialized (initialize the *b* coefficients) before used.
private
Expand Down Expand Up @@ -100,22 +101,157 @@ elemental subroutine init(self, steps)
select case(steps)
case(1)
! AB(1) Forward-Euler
self%b(1) = 1._R_P
self%b(1) = 1.0_R_P
case(2)
! AB(2)
self%b(1) = -0.5_R_P
self%b(2) = 1.5_R_P
self%b(1) = -1.0_R_P/2.0_R_P
self%b(2) = 3.0_R_P/2.0_R_P
case(3)
! AB(3)
self%b(1) = 5._R_P/12._R_P
self%b(2) = -4._R_P/3._R_P
self%b(3) = 23._R_P/12._R_P
self%b(1) = 5.0_R_P/12.0_R_P
self%b(2) = -16.0_R_P/12.0_R_P
self%b(3) = 23.0_R_P/12.0_R_P
case(4)
! AB(4)
self%b(1) = -3._R_P/8._R_P
self%b(2) = 37._R_P/24._R_P
self%b(3) = -59._R_P/24._R_P
self%b(4) = 55._R_P/24._R_P
self%b(1) = -9.0_R_P/24.0_R_P
self%b(2) = 37.0_R_P/24.0_R_P
self%b(3) = -59.0_R_P/24.0_R_P
self%b(4) = 55.0_R_P/24.0_R_P
case(5)
self%b(1) = 251.0_R_P/720.0_R_P
self%b(2) = -1274.0_R_P/720.0_R_P
self%b(3) = 2616.0_R_P/720.0_R_P
self%b(4) = -2774.0_R_P/720.0_R_P
self%b(5) = 1901.0_R_P/720.0_R_P
case(6)
self%b(1) = -475.0_R_P/1440.0_R_P
self%b(2) = 2877.0_R_P/1440.0_R_P
self%b(3) = -7298.0_R_P/1440.0_R_P
self%b(4) = 9982.0_R_P/1440.0_R_P
self%b(5) = -7923.0_R_P/1440.0_R_P
self%b(6) = 4277.0_R_P/1440.0_R_P
case(7)
self%b(1) = 19087.0_R_P/60480.0_R_P
self%b(2) = -134472.0_R_P/60480.0_R_P
self%b(3) = 407139.0_R_P/60480.0_R_P
self%b(4) = -688256.0_R_P/60480.0_R_P
self%b(5) = 705549.0_R_P/60480.0_R_P
self%b(6) = -447288.0_R_P/60480.0_R_P
self%b(7) = 198721.0_R_P/60480.0_R_P
case(8)
self%b(1) = -36799.0_R_P/120960.0_R_P
self%b(2) = 295767.0_R_P/120960.0_R_P
self%b(3) = -1041723.0_R_P/120960.0_R_P
self%b(4) = 2102243.0_R_P/120960.0_R_P
self%b(5) = -2664477.0_R_P/120960.0_R_P
self%b(6) = 2183877.0_R_P/120960.0_R_P
self%b(7) = -1152169.0_R_P/120960.0_R_P
self%b(8) = 434241.0_R_P/120960.0_R_P
case(9)
self%b(1) = 1070017.0_R_P/3628800.0_R_P
self%b(2) = -9664106.0_R_P/3628800.0_R_P
self%b(3) = 38833486.0_R_P/3628800.0_R_P
self%b(4) = -91172642.0_R_P/3628800.0_R_P
self%b(5) = 137968480.0_R_P/3628800.0_R_P
self%b(6) = -139855262.0_R_P/3628800.0_R_P
self%b(7) = 95476786.0_R_P/3628800.0_R_P
self%b(8) = -43125206.0_R_P/3628800.0_R_P
self%b(9) = 14097247.0_R_P/3628800.0_R_P
case(10)
self%b(1) = -2082753.0_R_P/7257600.0_R_P
self%b(2) = 20884811.0_R_P/7257600.0_R_P
self%b(3) = -94307320.0_R_P/7257600.0_R_P
self%b(4) = 252618224.0_R_P/7257600.0_R_P
self%b(5) = -444772162.0_R_P/7257600.0_R_P
self%b(6) = 538363838.0_R_P/7257600.0_R_P
self%b(7) = -454661776.0_R_P/7257600.0_R_P
self%b(8) = 265932680.0_R_P/7257600.0_R_P
self%b(9) = -104995189.0_R_P/7257600.0_R_P
self%b(10) = 30277247.0_R_P/7257600.0_R_P
case(11)
self%b(1) = 134211265.0_R_P/479001600.0_R_P
self%b(2) = -1479574348.0_R_P/479001600.0_R_P
self%b(3) = 7417904451.0_R_P/479001600.0_R_P
self%b(4) = -22329634920.0_R_P/479001600.0_R_P
self%b(5) = 44857168434.0_R_P/479001600.0_R_P
self%b(6) = -63176201472.0_R_P/479001600.0_R_P
self%b(7) = 63716378958.0_R_P/479001600.0_R_P
self%b(8) = -46113029016.0_R_P/479001600.0_R_P
self%b(9) = 23591063805.0_R_P/479001600.0_R_P
self%b(10) = -8271795124.0_R_P/479001600.0_R_P
self%b(11) = 2132509567.0_R_P/479001600.0_R_P
case(12)
self%b(1) = -262747265.0_R_P/958003200.0_R_P
self%b(2) = 3158642445.0_R_P/958003200.0_R_P
self%b(3) = -17410248271.0_R_P/958003200.0_R_P
self%b(4) = 58189107627.0_R_P/958003200.0_R_P
self%b(5) = -131365867290.0_R_P/958003200.0_R_P
self%b(6) = 211103573298.0_R_P/958003200.0_R_P
self%b(7) = -247741639374.0_R_P/958003200.0_R_P
self%b(8) = 214139355366.0_R_P/958003200.0_R_P
self%b(9) = -135579356757.0_R_P/958003200.0_R_P
self%b(10) = 61633227185.0_R_P/958003200.0_R_P
self%b(11) = -19433810163.0_R_P/958003200.0_R_P
self%b(12) = 4527766399.0_R_P/958003200.0_R_P
case(13)
self%b(1) = 703604254357.0_R_P/2615348736000.0_R_P
self%b(2) = -9160551085734.0_R_P/2615348736000.0_R_P
self%b(3) = 55060974662412.0_R_P/2615348736000.0_R_P
self%b(4) = -202322913738370.0_R_P/2615348736000.0_R_P
self%b(5) = 507140369728425.0_R_P/2615348736000.0_R_P
self%b(6) = -915883387152444.0_R_P/2615348736000.0_R_P
self%b(7) = 1226443086129408.0_R_P/2615348736000.0_R_P
self%b(8) = -1233589244941764.0_R_P/2615348736000.0_R_P
self%b(9) = 932884546055895.0_R_P/2615348736000.0_R_P
self%b(10) = -524924579905150.0_R_P/2615348736000.0_R_P
self%b(11) = 214696591002612.0_R_P/2615348736000.0_R_P
self%b(12) = -61497552797274.0_R_P/2615348736000.0_R_P
self%b(13) = 13064406523627.0_R_P/2615348736000.0_R_P
case(14)
self%b(1) = -1382741929621.0_R_P/5230697472000.0_R_P
self%b(2) = 19382853593787.0_R_P/5230697472000.0_R_P
self%b(3) = -126174972681906.0_R_P/5230697472000.0_R_P
self%b(4) = 505586141196430.0_R_P/5230697472000.0_R_P
self%b(5) = -1393306307155755.0_R_P/5230697472000.0_R_P
self%b(6) = 2793869602879077.0_R_P/5230697472000.0_R_P
self%b(7) = -4204551925534524.0_R_P/5230697472000.0_R_P
self%b(8) = 4825671323488452.0_R_P/5230697472000.0_R_P
self%b(9) = -4246767353305755.0_R_P/5230697472000.0_R_P
self%b(10) = 2854429571790805.0_R_P/5230697472000.0_R_P
self%b(11) = -1445313351681906.0_R_P/5230697472000.0_R_P
self%b(12) = 537247052515662.0_R_P/5230697472000.0_R_P
self%b(13) = -140970750679621.0_R_P/5230697472000.0_R_P
self%b(14) = 27511554976875.0_R_P/5230697472000.0_R_P
case(15)
self%b(1) = 8164168737599.0_R_P/31384184832000.0_R_P
self%b(2) = -122594813904112.0_R_P/31384184832000.0_R_P
self%b(3) = 859236476684231.0_R_P/31384184832000.0_R_P
self%b(4) = -3728807256577472.0_R_P/31384184832000.0_R_P
self%b(5) = 11205849753515179.0_R_P/31384184832000.0_R_P
self%b(6) = -24704503655607728.0_R_P/31384184832000.0_R_P
self%b(7) = 41280216336284259.0_R_P/31384184832000.0_R_P
self%b(8) = -53246738660646912.0_R_P/31384184832000.0_R_P
self%b(9) = 53471026659940509.0_R_P/31384184832000.0_R_P
self%b(10) = -41825269932507728.0_R_P/31384184832000.0_R_P
self%b(11) = 25298910337081429.0_R_P/31384184832000.0_R_P
self%b(12) = -11643637530577472.0_R_P/31384184832000.0_R_P
self%b(13) = 3966421670215481.0_R_P/31384184832000.0_R_P
self%b(14) = -960122866404112.0_R_P/31384184832000.0_R_P
self%b(15) = 173233498598849.0_R_P/31384184832000.0_R_P
case(16)
self%b(1) = -16088129229375.0_R_P/62768369664000.0_R_P
self%b(2) = 257650275915823.0_R_P/62768369664000.0_R_P
self%b(3) = -1934443196892599.0_R_P/62768369664000.0_R_P
self%b(4) = 9038571752734087.0_R_P/62768369664000.0_R_P
self%b(5) = -29417910911251819.0_R_P/62768369664000.0_R_P
self%b(6) = 70724351582843483.0_R_P/62768369664000.0_R_P
self%b(7) = -129930094104237331.0_R_P/62768369664000.0_R_P
self%b(8) = 186087544263596643.0_R_P/62768369664000.0_R_P
self%b(9) = -210020588912321949.0_R_P/62768369664000.0_R_P
self%b(10) = 187463140112902893.0_R_P/62768369664000.0_R_P
self%b(11) = -131963191940828581.0_R_P/62768369664000.0_R_P
self%b(12) = 72558117072259733.0_R_P/62768369664000.0_R_P
self%b(13) = -30607373860520569.0_R_P/62768369664000.0_R_P
self%b(14) = 9622096909515337.0_R_P/62768369664000.0_R_P
self%b(15) = -2161567671248849.0_R_P/62768369664000.0_R_P
self%b(16) = 362555126427073.0_R_P/62768369664000.0_R_P
endselect
return
!---------------------------------------------------------------------------------------------------------------------------------
Expand Down
6 changes: 3 additions & 3 deletions src/lib/foodie_integrator_adams_bashforth_moulton.f90
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
!< FOODIE integrator: provide a predictor-corrector class of Adams-Bashforth-Moutlon multi-step schemes, from 1st to 4rd order
!< FOODIE integrator: provide a predictor-corrector class of Adams-Bashforth-Moutlon multi-step schemes, from 1st to 16th order
!< accurate.
module foodie_integrator_adams_bashforth_moulton
!-----------------------------------------------------------------------------------------------------------------------------------
!< FOODIE integrator: provide a predictor-corrector class of Adams-Bashforth-Moutlon multi-step schemes, from 1st to 4rd order
!< FOODIE integrator: provide a predictor-corrector class of Adams-Bashforth-Moutlon multi-step schemes, from 1st to 16th order
!< accurate.
!<
!< Considering the following ODE system:
Expand Down Expand Up @@ -92,7 +92,7 @@ module foodie_integrator_adams_bashforth_moulton

!-----------------------------------------------------------------------------------------------------------------------------------
type :: adams_bashforth_moulton_integrator
!< FOODIE integrator: provide an explicit class of Adams-Bashforth-Moulton multi-step schemes, from 1st to 4rd order accurate.
!< FOODIE integrator: provide an explicit class of Adams-Bashforth-Moulton multi-step schemes, from 1st to 16th order accurate.
!<
!< @note The integrator must be created or initialized (predictor and corrector schemes selection) before used.
private
Expand Down

2 comments on commit 360c25a

@szaghi
Copy link

@szaghi szaghi commented on 360c25a Dec 30, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wonderful! @jacobwilliams very great job! As soon as possible I will grab these new coefficients, but if you like to make a pull request from your fork it is very welcome .

I never used/studied multi-step methods of such an order: the stability locus of such schemes coukd be an issue... does how change their stability?

@jacobwilliams
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll do a pull request after I test it some more. Yes, they can be unstable at the very high orders (actually this can be seen in your oscillation test case). For double precision, I've used up to 8th order for fixed step ABM (and also variable step ddeabm goes up to 12th order).

Please sign in to comment.