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