-
Notifications
You must be signed in to change notification settings - Fork 1
/
solvers.f90
239 lines (202 loc) · 13.8 KB
/
solvers.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
module solvers
use precision
use functions
implicit none
private
public :: rk2
public :: rk4
public :: dopri54
public :: dopri87
public :: solver
interface
subroutine solver(ff, t, dt, u0, u, err)
use precision
interface
function ff(t,u) result(up)
use precision
real(dp), intent(in) :: t
real(dp), intent(in) :: u(:)
real(dp), allocatable :: up(:)
end function
end interface
real(dp), intent(in) :: t
real(dp), intent(in) :: dt
real(dp), intent(in) :: u0(:)
real(dp), intent(inout) :: u(:)
real(dp), intent(inout) :: err
end subroutine solver
end interface
contains
subroutine rk2(f, t, dt, u0, u)
procedure(func) :: f
real(dp), intent(in) :: t
real(dp), intent(in) :: dt
real(dp), intent(in) :: u0(:)
real(dp), intent(inout) :: u(:)
real(dp), allocatable :: k1(:), k2(:)
k1 = f(t,u0)
k2 = f(t+dt, u0 + dt*k1)
u = u0 + (k1+k2)*dt*0.5_dp
end subroutine rk2
subroutine rk4(f, t, dt, u0, u)
procedure(func) :: f
real(dp), intent(in) :: t
real(dp), intent(in) :: dt
real(dp), intent(in) :: u0(:)
real(dp), intent(inout) :: u(:)
real(dp), allocatable :: k1(:), k2(:), k3(:), k4(:)
k1 = f(t , u0 )
k2 = f(t+dt*0.5_dp, u0+dt*k1*0.5_dp)
k3 = f(t+dt*0.5_dp, u0+dt*k2*0.5_dp)
k4 = f(t+dt , u0+dt*k3 )
u = u0 + (0.16666666666666666666_dp*k1+0.33333333333333333333_dp*k2+ &
& 0.33333333333333333333_dp*k3+0.16666666666666666666_dp*k4)*dt
end subroutine rk4
! --------------------------------------------------
subroutine dopri54(f, t, h, u0, u, err)
procedure(func) :: f
real(dp), intent(in) :: t, h
real(dp), dimension(:), intent(in) :: u0
real(dp), dimension(:), intent(inout) :: u
real(dp), intent(inout) :: err
real(dp), dimension(:,:), allocatable :: k
allocate(k(size(u0) ,7))
k(:,1) = f(t, u0)*h
k(:,2) = f(t+h/5.0_dp, u0+k(:,1)/5.0_dp)*h
k(:,3) = f(t+3.0_dp*h/10.0_dp, u0+3.0_dp/40.0_dp*k(:,1)+9.0_dp/40.0_dp*k(:,2))*h
k(:,4) = f(t+h*4.0_dp/5.0_dp, u0+44.0_dp/45.0_dp*k(:,1)+160.0_dp/45.0_dp*k(:,3)-168.0_dp/45.0_dp*k(:,2))*h
k(:,5) = f(t+h*8.d0/9.d0, u0+38744.d0/13122.d0*k(:,1)+128896.d0/13122.d0*k(:,3)-152160.d0/13122.d0*k(:,2) &
-3816.d0/13122.d0*k(:,4))*h
k(:,6) = f(t+h, u0+9017.d0/3168.d0*k(:,1)+46732.d0/5247.d0*k(:,3)+5194.d0/18656.d0*k(:,4) &
-34080.d0/3168.d0*k(:,2)-5103.d0/18656.d0*k(:,5))*h
u = u0+(35.d0/384.d0*k(:,1)+500.d0/1113.d0*k(:,3)+125.d0/192.d0*k(:,4)+&
11.d0/84.d0*k(:,6)-2187.d0/6784.d0*k(:,5))
k(:,7) = f(t+h, u)*h
! just estimate the error as ( z - y ):
! err = max(|u(5)-u(4)|)
! maxval(u(:)) maxloc(u(:)) indice del valore massimo
err = maxval(abs( 71.d0/57600.d0*k(:,1)+22.d0/525.d0*k(:,6) &
& +71.d0/1920.d0*k(:,4) &
& -71.d0/16695.d0*k(:,3)-17253.d0/339200.d0*k(:,5)-1.d0/40.d0*k(:,7) ))
deallocate(k)
end subroutine dopri54
! --------------------------------------------------
subroutine dopri87(f, t, dt, u0, u, err)
procedure(func) :: f
real(dp), intent(in) :: t
real(dp), intent(in) :: dt
real(dp), intent(in) :: u0(:)
real(dp), intent(inout) :: u(:)
real(dp), intent(inout) :: err
! Parameters taken from
real(dp), parameter :: c1= 0.0_dp
real(dp), parameter :: c2=.05555555555555555555555555555555555555555555555555555555555555555555555555555555555556_dp
real(dp), parameter :: c3=.08333333333333333333333333333333333333333333333333333333333333333333333333333333333333_dp
real(dp), parameter :: c4=.125_dp
real(dp), parameter :: c5=.3125_dp
real(dp), parameter :: c6=.375_dp
real(dp), parameter :: c7=.1475_dp
real(dp), parameter :: c8=.465_dp
real(dp), parameter :: c9=.5648654513822595753983585014261682587385670087264133215107858541861772043626893634869_dp
real(dp), parameter :: c10=.65_dp
real(dp), parameter :: c11=.9246562776405044467450135743183695426492034467027398177693057974193090932038137681734_dp
real(dp), parameter :: c12=1.0_dp
real(dp), parameter :: c13=1.0_dp
real(dp), parameter :: a21=.05555555555555555555555555555555555555555555555555555555555555555555555555555555555556_dp
real(dp), parameter :: a31=.02083333333333333333333333333333333333333333333333333333333333333333333333333333333333_dp
real(dp), parameter :: a32=.0625_dp
real(dp), parameter :: a41=.03125_dp
real(dp), parameter :: a43=.09375_dp
real(dp), parameter :: a51=.3125_dp
real(dp), parameter :: a53=-1.171875_dp
real(dp), parameter :: a54=1.171875_dp
real(dp), parameter :: a61=.0375_dp
real(dp), parameter :: a64=.1875_dp
real(dp), parameter :: a65=.15_dp
real(dp), parameter :: a71=.04791013711111111111111111111111111111111111111111111111111111111111111111111111111111_dp
real(dp), parameter :: a74=.1122487127777777777777777777777777777777777777777777777777777777777777777777777777778_dp
real(dp), parameter :: a75=-.02550567377777777777777777777777777777777777777777777777777777777777777777777777777778_dp
real(dp), parameter :: a76=.01284682388888888888888888888888888888888888888888888888888888888888888888888888888889_dp
real(dp), parameter :: a81=.01691798978729228118143110713603823606551492879543441068765183916901985069878116840258_dp
real(dp), parameter :: a84=.3878482784860431695265457441593733533707275558526027137301504136188413357699752956243_dp
real(dp), parameter :: a85=.03597736985150032789670088963477236800815873945874968299566918115320174598617642813621_dp
real(dp), parameter :: a86=.1969702142156660601567152560721498881281698021684239074332818272245101975612112344987_dp
real(dp), parameter :: a87=-.1727138523405018387613929970023338455725710262752107148467532611655731300161441266618_dp
real(dp), parameter :: a91=.06909575335919230064856454898454767856104137244685570409618590205329379141801779261271_dp
real(dp), parameter :: a94=-.6342479767288541518828078749717385465426336031120747443029278238678259156558727343870_dp
real(dp), parameter :: a95=-.1611975752246040803668769239818171234422237864808493434053355718725564004275765826294_dp
real(dp), parameter :: a96=.1386503094588252554198669501330158019276654889495806914244800957316809692178564181463_dp
real(dp), parameter :: a97=.9409286140357562697242396841302583434711811483145028697758469500622977999086075976730_dp
real(dp), parameter :: a98=.2116363264819439818553721171319021047635363886083981439225363020792869599016568720713_dp
real(dp), parameter :: a101=.1835569968390453854898060235368803084497642516277329033567822939550366893470341427061_dp
real(dp), parameter :: a104=-2.468768084315592452744315759974107457777649753547732495587510208118948846864075671103_dp
real(dp), parameter :: a105=-.2912868878163004563880025728039519800543380294081727678722180378521228988619909525814_dp
real(dp), parameter :: a106=-.02647302023311737568843979946594614325963055282676866851254669985184857900430792761417_dp
real(dp), parameter :: a107=2.847838764192800449164518254216773770231582185011953158945518598604677368771006006947_dp
real(dp), parameter :: a108=.2813873314698497925394036418267117820980705455360140173438806414700945017562775967023_dp
real(dp), parameter :: a109=.1237448998633146576270302126636397203122013536069738523260934117931117648560568049432_dp
real(dp), parameter :: a111=-1.215424817395888059160510525029662994880024988044112261492933435562347094075572196306_dp
real(dp), parameter :: a114=16.67260866594577243228041328856410774858907460078758981907659463146817850515953150318_dp
real(dp), parameter :: a115=.9157418284168179605957186504507426331593736334298114556675054711691705693180672015549_dp
real(dp), parameter :: a116=-6.056605804357470947554505543091634004081967083324413691952297768849489422976536037188_dp
real(dp), parameter :: a117=-16.00357359415617811184170641007882303068079304063907122528682690677051264091851070681_dp
real(dp), parameter :: a118=14.84930308629766255754539189802663208272299893302745636963476380528935538206408294404_dp
real(dp), parameter :: a119=-13.37157573528984931829304139618159579089195289469740214314819172008509426917249413996_dp
real(dp), parameter :: a1110=5.134182648179637933173253611658602898712494286162881495270691720760048063805245199662_dp
real(dp), parameter :: a121=.2588609164382642838157309322317577667296307766301063163257807024364127266506000943372_dp
real(dp), parameter :: a124=-4.774485785489205112310117509706042746829391853746564335432052648850273207666153414325_dp
real(dp), parameter :: a125=-.4350930137770325094407004118103177819323551661617974116300885410959889587870529265443_dp
real(dp), parameter :: a126=-3.049483332072241509560512866312031613982854911220735245095857885009959537455491093975_dp
real(dp), parameter :: a127=5.577920039936099117423676634464941858623588944531498679305832747426169795721583100328_dp
real(dp), parameter :: a128=6.155831589861040097338689126688954481197754937462937466615262374558053507207577588692_dp
real(dp), parameter :: a129=-5.062104586736938370077406433910391644990220712141673880668482097753119682588189892664_dp
real(dp), parameter :: a1210=2.193926173180679061274914290465806019788262707389033759251111926114017568440058142981_dp
real(dp), parameter :: a1211=.1346279986593349415357262378873236613955852772571946513284934221746877884770684011715_dp
real(dp), parameter :: a131=.8224275996265074779631682047726665909572303617765850630130165537026064642659285212794_dp
real(dp), parameter :: a134=-11.65867325727766428397655303545841477547369082638864247596731917545919361630644756876_dp
real(dp), parameter :: a135=-.7576221166909361958811161540882449653663757591954118634754438929149012424414834787423_dp
real(dp), parameter :: a136=.7139735881595815279782692827650546753142248878566301594716125352788068216039494192152_dp
real(dp), parameter :: a137=12.07577498689005673956617044860067967095705800972197897084832370633043082304810801069_dp
real(dp), parameter :: a138=-2.127659113920402656390820858969398635427927973275119750336406473203376323696576615278_dp
real(dp), parameter :: a139=1.990166207048955418328071698344314152176173017359794982793519250552799420955298804232_dp
real(dp), parameter :: a1310=-.2342864715440402926602946918568015314512425817004394700255034276765120670064339827205_dp
real(dp), parameter :: a1311=.1758985777079422650731051058901448183145508638446243836782009233893397195776568900819_dp
real(dp), parameter :: b81=.04174749114153024622208592846850711513419360280744778485846131359556941056165592859127_dp
real(dp), parameter :: b86=-.05545232861123930896152189465471671889359328156519226303661201143086248409346173220672_dp
real(dp), parameter :: b87=.2393128072011800970467473542487569696603053593110659071772286813575015789696127063297_dp
real(dp), parameter :: b88=.7035106694034430230580464108897021513663799729422329988912282232809616439130702896567_dp
real(dp), parameter :: b89=-.7597596138144609298844876770850584076554230192157088121727371293320966312254783676226_dp
real(dp), parameter :: b810=.6605630309222863414613785948378206399404197125341442757648140724130090459250911959182_dp
real(dp), parameter :: b811=.1581874825101233355296148386006854439728567324034642678211136427219891617064005758134_dp
real(dp), parameter :: b812=-.2381095387528628044718635553056971935251390792174541593034967926060717257568905964800_dp
real(dp), parameter :: b813=.25_dp
real(dp), parameter :: b71=.02955321367635349698196488311203224657733277925009198076731622388624965866492863628737_dp
real(dp), parameter :: b76=-.8286062764877970397668056126887191847354037544018122961965538027330129339209559093501_dp
real(dp), parameter :: b77=.3112409000511183279299137516268570512897363217826855582698888340873138124416908811814_dp
real(dp), parameter :: b78=2.467345190599886981964685704068761458562259575475834628418174165533324998199993514039_dp
real(dp), parameter :: b79=-2.546941651841908739127380075415708961780905163114681004322493329883513805107784462867_dp
real(dp), parameter :: b710=1.443548583676775240301874950690104268511068010189240811834537552274685797398657155566_dp
real(dp), parameter :: b711=.07941559588112728727130195416222867713146778637419587678468591239050802787902574069868_dp
real(dp), parameter :: b712=.04444444444444444444444444444444444444444444444444444444444444444444444444444444444444_dp
real(dp), allocatable :: k1(:), k2(:), k3(:), k4(:), k5(:), k6(:)
real(dp), allocatable :: k7(:), k8(:), k9(:), k10(:), k11(:), k12(:), k13(:)
real(dp), allocatable :: y(:)
k1 = dt*f(t , u0 )
k2 = dt*f(t+dt*c2 , u0+ a21*k1)
k3 = dt*f(t+dt*c3 , u0+ a31*k1+ a32*k2)
k4 = dt*f(t+dt*c4 , u0+ a41*k1+ a43*k3)
k5 = dt*f(t+dt*c5 , u0+ a51*k1+ a53*k3+ a54*k4)
k6 = dt*f(t+dt*c6 , u0+ a61*k1+ a64*k4+ a65*k5)
k7 = dt*f(t+dt*c7 , u0+ a71*k1+ a74*k4+ a75*k5+ a76*k6)
k8 = dt*f(t+dt*c8 , u0+ a81*k1+ a84*k4+ a85*k5+ a86*k6+ a87*k7)
k9 = dt*f(t+dt*c9 , u0+ a91*k1+ a94*k4+ a95*k5+ a96*k6+ a97*k7+ a98*k8)
k10= dt*f(t+dt*c10, u0+a101*k1+a104*k4+a105*k5+a106*k6+a107*k7+a108*k8+a109*k9)
k11= dt*f(t+dt*c11, u0+a111*k1+a114*k4+a115*k5+a116*k6+a117*k7+a118*k8+a119*k9+a1110*k10)
k12= dt*f(t+dt*c12, u0+a121*k1+a124*k4+a125*k5+a126*k6+a127*k7+a128*k8+a129*k9+a1210*k10+a1211*k11)
k13= dt*f(t+dt*c13, u0+a131*k1+a134*k4+a135*k5+a136*k6+a137*k7+a138*k8+a139*k9+a1310*k10+a1311*k11)
u = u0 + b81*k1+b86*k6+b87*k7+b88*k8+b89*k9+b810*k10+b811*k11+b812*k12+b813*k13
! Error estimation |u-y|
y = u0 + b71*k1+b76*k6+b77*k7+b78*k8+b79*k9+b710*k10+b711*k11+b712*k12
err = maxval(abs(u-y))
end subroutine dopri87
end module solvers