Skip to content

Commit

Permalink
added rktf65
Browse files Browse the repository at this point in the history 
working on another test case
  • Loading branch information
@jacobwilliams
jacobwilliams committed Jul 5, 2023
1 parent 7032b45 commit 817f9e4
Show file tree
Hide file tree
Showing 5 changed files with 352 additions and 7 deletions.
3 changes: 2 additions & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ Method name | Type | Order | Number of Stages | Reference
`rktp64` | Variable-step | 6 | 7 | [Tsitouras & Papakostas (1999)](https://epubs.siam.org/doi/abs/10.1137/S1064827596302230?journalCode=sjoce3)
`rkv65e` | Variable-step | 6 | 9 (FSAL) | [Verner (1994)](https://www.sfu.ca/~jverner/RKV65.IIIXb.Efficient.00000144617.081204.CoeffsOnlyFLOAT)
`rkv65r` | Variable-step | 6 | 9 (FSAL) | [Verner (1994)](https://www.sfu.ca/~jverner/RKV65.IIIXb.Robust.00010102836.081204.RATOnWeb)
`rktf65` | Variable-step | 6 | 9 (FSAL) | [Tsitouras & Famelis (2006)](http://users.uoa.gr/~tsitourasc/ModifiedRK-ICNAAM2006.pdf)
`rktp75` | Variable-step | 7 | 9 | [Tsitouras & Papakostas (1999)](https://epubs.siam.org/doi/abs/10.1137/S1064827596302230?journalCode=sjoce3)
`rktmy7` | Variable-step | 7 | 10 | [Tanaka, Muramatsu & Yamashita (1992)](http://www.peterstone.name/Maplepgs/Maple/nmthds/RKcoeff/Runge_Kutta_schemes/RK7/RKcoeff7d_4.pdf)
`rkv76e` | Variable-step | 7 | 10 | [Verner (1978)](https://epubs.siam.org/doi/10.1137/0715051)
Expand All @@ -67,7 +68,7 @@ Method name | Type | Order | Number of Stages | Reference
`rkv98r` | Variable-step | 9 | 16 | [Verner (1978)](https://www.jstor.org/stable/2156853)
`rkf108` | Variable-step | 10 | 17 | [Feagin (2006)](https://sce.uhcl.edu/rungekutta/rk108.txt)
`rkc108` | Variable-step | 10 | 21 | [Curtis (1975)](http://www.peterstone.name/Maplepgs/Maple/nmthds/RKcoeff/Runge_Kutta_schemes/RK10/RKcoeff10a(8)_2.pdf)
`rks1110a` | Variable-step | 11 | 26 | [Stone (2015)](http://www.peterstone.name/Maplepgs/Maple/nmthds/RKcoeff/Runge_Kutta_schemes/RK11/RKcoeff11_a.pdf)
`rks1110a` | Variable-step | 11 | 26 | [Stone (2015)](http://www.peterstone.name/Maplepgs/Maple/nmthds/RKcoeff/Runge_Kutta_schemes/RK11/RKcoeff11_a.pdf)
`rkf1210` | Variable-step | 12 | 25 | [Feagin (2006)](https://sce.uhcl.edu/rungekutta/rk1210.txt)
`rko129` | Variable-step | 12 | 29 | [Ono (2006)](http://www.peterstone.name/Maplepgs/Maple/nmthds/RKcoeff/Runge_Kutta_schemes/RK12/RKcoeff12h(9)_1.pdf)
`rkf1412` | Variable-step | 14 | 35 | [Feagin (2006)](https://sce.uhcl.edu/rungekutta/rk1412.txt)
Expand Down
32 changes: 26 additions & 6 deletions src/rklib_module.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -270,25 +270,31 @@ module rklib_module
procedure :: order => rktp64_order
end type rktp64_class
type,extends(rk_variable_step_fsal_class),public :: rkv65e_class
!! Verner's 'most efficient' Runge--Kutta (9,6(5))
!! Verner's 'most efficient' Runge-Kutta (9,6(5))
contains
procedure :: step => rkv65e
procedure :: order => rkv65e_order
end type rkv65e_class
type,extends(rk_variable_step_fsal_class),public :: rktf65_class
!! Tsitouras & Famelis 6(5)
contains
procedure :: step => rktf65
procedure :: order => rktf65_order
end type rktf65_class
type,extends(rk_variable_step_fsal_class),public :: rkv65r_class
!! Verner's 'most robust' Runge--Kutta (9,6(5))
!! Verner's 'most robust' Runge-Kutta (9,6(5))
contains
procedure :: step => rkv65r
procedure :: order => rkv65r_order
end type rkv65r_class
type,extends(rk_variable_step_class),public :: rkv76e_class
!! Verner's 'most efficient' Runge--Kutta (10:7(6))
!! Verner's 'most efficient' Runge-Kutta (10:7(6))
contains
procedure :: step => rkv76e
procedure :: order => rkv76e_order
end type rkv76e_class
type,extends(rk_variable_step_class),public :: rkv76r_class
!! Verner's 'most robust' Runge--Kutta (10:7(6))
!! Verner's 'most robust' Runge-Kutta (10:7(6))
contains
procedure :: step => rkv76r
procedure :: order => rkv76r_order
Expand All @@ -306,13 +312,13 @@ module rklib_module
procedure :: order => rkdp87_order
end type rkdp87_class
type,extends(rk_variable_step_class),public :: rkv87e_class
!! Verner's "most efficient" Runge--Kutta (8)7 method.
!! Verner's "most efficient" Runge-Kutta (8)7 method.
contains
procedure :: step => rkv87e
procedure :: order => rkv87e_order
end type rkv87e_class
type,extends(rk_variable_step_class),public :: rkv87r_class
!! Verner's "most robust" Runge--Kutta (8)7 method.
!! Verner's "most robust" Runge-Kutta (8)7 method.
contains
procedure :: step => rkv87r
procedure :: order => rkv87r_order
Expand Down Expand Up @@ -688,6 +694,15 @@ module subroutine rkv65e(me,t,x,h,xf,terr)
real(wp),dimension(me%n),intent(out) :: xf
real(wp),dimension(me%n),intent(out) :: terr
end subroutine rkv65e
module subroutine rktf65(me,t,x,h,xf,terr)
implicit none
class(rktf65_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(me%n),intent(in) :: x
real(wp),intent(in) :: h
real(wp),dimension(me%n),intent(out) :: xf
real(wp),dimension(me%n),intent(out) :: terr
end subroutine rktf65
module subroutine rkv65r(me,t,x,h,xf,terr)
implicit none
class(rkv65r_class),intent(inout) :: me
Expand Down Expand Up @@ -931,6 +946,11 @@ pure module function rkv65e_order(me) result(p)
class(rkv65e_class),intent(in) :: me
integer :: p !! order of the method
end function rkv65e_order
pure module function rktf65_order(me) result(p)
implicit none
class(rktf65_class),intent(in) :: me
integer :: p !! order of the method
end function rktf65_order
pure module function rkv65r_order(me) result(p)
implicit none
class(rkv65r_class),intent(in) :: me
Expand Down
104 changes: 104 additions & 0 deletions src/rklib_variable_submodule.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -684,6 +684,109 @@
end procedure rkv65e
!*****************************************************************************************

!*****************************************************************************************
!>
! Tsitouras & Famelis Runge-Kutta 6(5) method.
!
!### Reference
! * Ch. Tsitouras and I. Th. Famelis,
! [Phase-Fitted modified Runge-Kutta pairs of orders 6(5)](https://www.researchgate.net/publication/251740152_Phase-Fitted_modified_Runge-Kutta_pairs_of_orders_65),
! ICNAAM 2006, Crete, 2006.
! * [More precise rational coefficients](http://www.peterstone.name/Maplepgs/Maple/nmthds/RKcoeff/Runge_Kutta_schemes/RK6/RKcoeff6g_5.pdf)
!
!@note The floating point coefficients here were generated from the rational ones from the 2nd reference.
!@note This is a first-same-as-last (FSAL) step.

module procedure rktf65

real(wp),parameter :: a2 = 9.28961748633879781420765027322404371584699453551912568306010928961748633879781420765e-2_wp
real(wp),parameter :: a3 = 1.44578313253012048192771084337349397590361445783132530120481927710843373493975903614e-1_wp
real(wp),parameter :: a4 = 2.16867469879518072289156626506024096385542168674698795180722891566265060240963855422e-1_wp
real(wp),parameter :: a5 = 5.68000000000000000000000000000000000000000000000000000000000000000000000000000000000e-1_wp
real(wp),parameter :: a6 = 7.11864406779661016949152542372881355932203389830508474576271186440677966101694915254e-1_wp
real(wp),parameter :: a7 = 9.95000000000000000000000000000000000000000000000000000000000000000000000000000000000e-1_wp

real(wp),parameter :: b21 = 9.28961748633879781420765027322404371584699453551912568306010928961748633879781420765e-2_wp
real(wp),parameter :: b31 = 3.20715889781663863106572284887245651635599805316233039884555941697335052470690700435e-2_wp
real(wp),parameter :: b32 = 1.12506724274845661882113855848624832426801465251509226132026333541109868246906833571e-1_wp
real(wp),parameter :: b41 = 5.42168674698795180722891566265060240963855421686746987951807228915662650602409638554e-2_wp
real(wp),parameter :: b43 = 1.62650602409638554216867469879518072289156626506024096385542168674698795180722891566e-1_wp
real(wp),parameter :: b51 = 6.56598126617283950617283950617283950617283950617283950617283950617283950617283950617e-1_wp
real(wp),parameter :: b53 = -2.4972770465185185185185185185185185185185185185185185185185185185185185185185185185e0_wp
real(wp),parameter :: b54 = 2.40867891990123456790123456790123456790123456790123456790123456790123456790123456790e0_wp
real(wp),parameter :: b61 = -1.7212338830272756350112289300338046256920059840520750991645177506321959252208613366e0_wp
real(wp),parameter :: b63 = 7.22310751192021126636920155870658124799459070542113115974697356535677039490201029989e0_wp
real(wp),parameter :: b64 = -5.4959191727393518000813340726320306436524484390690989489777480799761402436788850115e0_wp
real(wp),parameter :: b65 = 7.05909950626077185672513986332135377282067107530551362971563451692243740099430963481e-1_wp
real(wp),parameter :: b71 = 4.12867189716161527469018071566157538750756968750429124243842122794843113666860197270e0_wp
real(wp),parameter :: b73 = -1.6914025304289438220281346004343739703895247887189472122417356767305063273262951589e1_wp
real(wp),parameter :: b74 = 1.43289813217405868844323023574174855877382144987680899683823277400972814702596244738e1_wp
real(wp),parameter :: b75 = -1.5533550903735309941352947451105003919180320410408681830389632181656533578398372284e0_wp
real(wp),parameter :: b76 = 1.00472717576076705529415767637517912056749574195795909463557101742500402417456237099e0_wp
real(wp),parameter :: b81 = 4.46860784284429599022330889710676757001193443328341367781153268526799285074027050787e0_wp
real(wp),parameter :: b83 = -1.8345441869429650066499014099391358577583747335752433112108627235199546393109090809e1_wp
real(wp),parameter :: b84 = 1.55237707233543199854527260736987323920859616527670939779256214402670438227794675418e1_wp
real(wp),parameter :: b85 = -1.7228800213316940986071454499455146268823602061289992205388429896934101951354882736e0_wp
real(wp),parameter :: b86 = 1.08151571740214517180934640319081789154027660669157008569533046194596896725782768317e0_wp
real(wp),parameter :: b87 = -5.5723928394169823792218246594446491720651508606454087850143625880490525329866504633e-3_wp
! real(wp),parameter :: b91 = 6.42309093721083246173486105540744934898802830008152394733186013531402447815498986849e-2_wp ! FSAL
! real(wp),parameter :: b94 = 3.32861824699421109595374929094856793240150741775741341700388991831307262888094726546e-1_wp
! real(wp),parameter :: b95 = 2.67859229165778080498089104036336501876704241846808785499560813816788199727204331902e-1_wp
! real(wp),parameter :: b96 = 1.79863899670938711343662526258687438581919457377762276997401953795124106054682426881e-1_wp
! real(wp),parameter :: b97 = 1.51075784805762160255845148219838318965956604578826959670696217369891977986727762107e0_wp
! real(wp),parameter :: b98 = -1.3555737109658678286129266521423384168482207697893972403776325344952795933188090051e0_wp

real(wp),parameter :: c1 = 6.42309093721083246173486105540744934898802830008152394733186013531402447815498986849e-2_wp
real(wp),parameter :: c4 = 3.32861824699421109595374929094856793240150741775741341700388991831307262888094726546e-1_wp
real(wp),parameter :: c5 = 2.67859229165778080498089104036336501876704241846808785499560813816788199727204331902e-1_wp
real(wp),parameter :: c6 = 1.79863899670938711343662526258687438581919457377762276997401953795124106054682426881e-1_wp
real(wp),parameter :: c7 = 1.51075784805762160255845148219838318965956604578826959670696217369891977986727762107e0_wp
real(wp),parameter :: c8 = -1.3555737109658678286129266521423384168482207697893972403776325344952795933188090051e0_wp

real(wp),parameter :: d1 = 6.22980954171238844288765383876964441563622264193747564518688416336360055072506207144e-2_wp
real(wp),parameter :: d4 = 3.40203519635783669970356076069438053031579043631460538976870898652474922063170157204e-1_wp
real(wp),parameter :: d5 = 2.35997541364108988248916956975923713134336326558729708330181566604979639253608698633e-1_wp
real(wp),parameter :: d6 = 2.20640276396042514942275961685396046511990069015617783025687162561106239011510198458e-1_wp
real(wp),parameter :: d7 = 1.15044133952868915480448599961883762022583301353216360785466101429534020477553812807e0_wp
real(wp),parameter :: d8 = -1.0029141056750815457282448660706252103934340124906797279726028170808703439444111364e0_wp
real(wp),parameter :: d9 = -6.6666666666666666666666666666666666666666666666666666666666666666666666666666666667e-3_wp

real(wp),parameter :: e1 = c1 - d1
real(wp),parameter :: e4 = c4 - d4
real(wp),parameter :: e5 = c5 - d5
real(wp),parameter :: e6 = c6 - d6
real(wp),parameter :: e7 = c7 - d7
real(wp),parameter :: e8 = c8 - d8
real(wp),parameter :: e9 = - d9

real(wp),dimension(me%n) :: f1,f2,f3,f4,f5,f6,f7,f8,f9

if (h==zero) then
xf = x
terr = zero
return
end if

! check the cached function eval of the last step:
call me%check_fsal_cache(t,x,f1)

call me%f(t+a2*h,x+h*(b21*f1),f2)
call me%f(t+a3*h,x+h*(b31*f1+b32*f2),f3)
call me%f(t+a4*h,x+h*(b41*f1+b43*f3),f4)
call me%f(t+a5*h,x+h*(b51*f1+b53*f3+b54*f4),f5)
call me%f(t+a6*h,x+h*(b61*f1+b63*f3+b64*f4+b65*f5),f6)
call me%f(t+a7*h,x+h*(b71*f1+b73*f3+b74*f4+b75*f5+b76*f6),f7)
call me%f(t+h, x+h*(b81*f1+b83*f3+b84*f4+b85*f5+b86*f6+b87*f7),f8)

! last point is cached for the next step:
xf = x+h*(c1*f1+c4*f4+c5*f5+c6*f6+c7*f7+c8*f8)
call me%set_fsal_cache(t+h,xf,f9)

terr = h*(e1*f1+e4*f4+e5*f5+e6*f6+e7*f7+e8*f8+e9*f9)

end procedure rktf65
!*****************************************************************************************

!*****************************************************************************************
!>
! Verner's "most robust" Runge-Kutta (9,6(5)) pair.
Expand Down Expand Up @@ -5988,6 +6091,7 @@
module procedure rkc65_order ; p = 6 ; end procedure !! Returns the order of the [[rkc65]] method.
module procedure rkv65e_order ; p = 6 ; end procedure !! Returns the order of the [[rkv65e]] method.
module procedure rkv65r_order ; p = 6 ; end procedure !! Returns the order of the [[rkv65r]] method.
module procedure rktf65_order ; p = 6 ; end procedure !! Returns the order of the [[rktf65]] method.
module procedure rktp75_order ; p = 7 ; end procedure !! Returns the order of the [[rktp75]] method.
module procedure rktmy7_order ; p = 7 ; end procedure !! Returns the order of the [[rktmy7]] method.
module procedure rkv76e_order ; p = 7 ; end procedure !! Returns the order of the [[rkv76e]] method.
Expand Down
Loading

0 comments on commit 817f9e4

Please sign in to comment.