Newer
Older
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
C_Coef = [0, .5, .5, .5 + np.sqrt(21)/ 14, .5 + np.sqrt(21)/ 14, .5, .5 - np.sqrt(21)/ 14, .5 - np.sqrt(21)/ 14, .5, .5 + np.sqrt(21)/ 14, 1]
B_Coef = [.05, 0, 0, 0, 0, 0, 0, 49/180, 16/45, 49/180, .05]
a = np.zeros((12,12))
a[2,1]=1/2
a[3,1]=1/4
a[3,2]=1/4
a[4,1]=1/7
a[4,2]=-1/14-3/98*21**(1/2)
a[4,3]=3/7+5/49*21**(1/2)
a[5,1]=11/84+1/84*21**(1/2)
a[5,2]=0
a[5,3]=2/7+4/63*21**(1/2)
a[5,4]=1/12-1/252*21**(1/2)
a[6,1]=5/48+1/48*21**(1/2)
a[6,2]=0
a[6,3]=1/4+1/36*21**(1/2)
a[6,4]=-77/120+7/180*21**(1/2)
a[6,5]=63/80-7/80*21**(1/2)
a[7,1]=5/21-1/42*21**(1/2)
a[7,2]=0
a[7,3]=-48/35+92/315*21**(1/2)
a[7,4]=211/30-29/18*21**(1/2)
a[7,5]=-36/5+23/14*21**(1/2)
a[7,6]=9/5-13/35*21**(1/2)
a[8,1]=1/14
a[8,2]=0
a[8,3]=0
a[8,4]=0
a[8,5]=1/9-1/42*21**(1/2)
a[8,6]=13/63-1/21*21**(1/2)
a[8,7]=1/9
a[9,1]=1/32
a[9,2]=0
a[9,3]=0
a[9,4]=0
a[9,5]=91/576-7/192*21**(1/2)
a[9,6]=11/72
a[9,7]=-385/1152-25/384*21**(1/2)
a[9,8]=63/128+13/128*21**(1/2)
a[10,1]=1/14
a[10,2]=0
a[10,3]=0
a[10,4]=0
a[10,5]=1/9
a[10,6]=-733/2205-1/15*21**(1/2)
a[10,7]=515/504+37/168*21**(1/2)
a[10,8]=-51/56-11/56*21**(1/2)
a[10,9]=132/245+4/35*21**(1/2)
a[11,1]=0
a[11,2]=0
a[11,3]=0
a[11,4]=0
a[11,5]=-7/3+7/18*21**(1/2)
a[11,6]=-2/5+28/45*21**(1/2)
a[11,7]=-91/24-53/72*21**(1/2)
a[11,8]=301/72+53/72*21**(1/2)
a[11,9]=28/45-28/45*21**(1/2)
a[11,10]=49/18-7/18*21**(1/2)
A_Coef = a[1:, 1:]
def k_func(k_step, func, t, state, delT, A_Coef, B_Coef, C_Coef, *args, **kwargs):
"""RK integration recursive function.
Args:
k_step: the current kfunction sub-step.
func: The ODE function to integrate.
t: The time of the initial condition state.
state: The value of the initial condition at t.
delT: The time step to integrate across.
A_Coef: The first set of Butcher Tableau Coefficients.
B_Coef: The second set of Butcher Tableau Coefficients.
C_Coef: The third set of Butcher Tableau Coefficients.
*args: Supplemental arguments for the ODE function func.
**kwargs: Supplemental argumentsfor ODE function func.
Returns:
The 'k_step'th k value for the integration step.
"""
if k_step == 1: #call the derivative function when necessary
return func(t, state, *args, **kwargs)
else:
carry = 0 #otherwise build the k value as a sum of the previous k values and thier coefficients recursively
for i in range(1, k_step):
carry += A_Coef[k_step - 1, i-1] * k_func(i-1, func, t, state, delT, A_Coef, B_Coef, C_Coef, *args, **kwargs)
return func(t + C_Coef[k_step - 1] * delT, state + delT*carry, *args, **kwargs)
def rk_step(func, t, state, delT, *args, **kwargs):
"""Generalized Runge-Kutta integrator
Args:
func: The ODE function to integrate.
t: The time of the initial condition state.
state: The value of the initial condition at t.
delT: The time step to integrate across.
A_Coef: The first set of Butcher Tableau Coefficients.
B_Coef: The second set of Butcher Tableau Coefficients.
C_Coef: The third set of Butcher Tableau Coefficients.
*args: Supplemental arguments for the ODE function func.
**kwargs: Supplemental argumentsfor ODE function func.
Returns:
Updated value of state at time t + delT.
"""
sum = 0 # sum of k functions and their coefficients.
for k_step in range(1, A_Coef.shape[1] + 1): #For each step according to the Butcher Tableau
sum += B_Coef[k_step - 1] * k_func(k_step, func, t, state, delT, A_Coef, B_Coef, C_Coef, *args, **kwargs) #Add contribution of each integration step
return state + delT * sum