I'm trying to implement a function that would input a butcher's tableau for a Runge-Kutta method of order $s$, an ode and its jacobian, as well as the initial values and return the next step. If we have a butcher tableau, $$ \begin{array}{c|cc} \alpha & \beta \\ \hline & \gamma \\ \end{array} $$ we'd like to convert that into a diagonally implicit butcher tableau $$ y_{n+1} = y_n + \sum_{i=1}^s \gamma_i \kappa_i $$ with \begin{align} \kappa_1 &= f\left(t_0 + \alpha_1 h, y_0 + \beta_{1j}\kappa_j\right) \\ \kappa_2 &= f\left(t_0 + \alpha_2 h, y_0 + h\sum_{j=1}^2 \beta_{2j}\kappa_j\right)\\ &\vdots \\ \kappa_s &= f\left(t_0 + \alpha_s h, y_0 + h\sum_{j=1}^s \beta_{sj}\kappa_j\right). \end{align} Keep in mind for the i-th step we need to solve for $\kappa_i$, where we use Newton Method to solve for a vector $\delta_i$ s.t. $F(\delta_i) = 0$, $$ F(\delta_i) = \delta_i - f\left(t_0 + \alpha_i h, y_0 + \beta_{ii}\delta_i + h\sum_{j=1}^{i-1} \beta_{ij}\kappa_j\right) $$ where DF is our Jacobian matrix.
def dirk(f, Df, t0, y0, h, alpha, beta, gamma): #---# Parameters #-------------------------------------------# #f(t,y) : lambda function of the ODE y'(t) = f(t, y) #Df(t,y): lambda function for the Jacobian of f #t0 : initial time #y0 : initial coordinate #alpha : Butcher tableau vector of size s (weights of y step) #beta : Butcher tableau matrix of size s*s (weights of steps) #gamma : Butcher tableau vector of size s (weights of k) #------------------------------------------------------------# n = len(beta) m = len(f(t0, y0)) K = np.zeros(len(alpha)) if beta[0][0] != 0: K0,_ = newton(lambda delta: delta - f(t0 + alpha[0]*h, y0+h*beta[0][0]*delta), lambda delta: np.eye(len(y0))-Df(t0+alpha[0]*h,y0+h*beta[0][0]*delta), f(t0, y0), 1e-8,1000) else: K0 = y0 + h*f(t0, y0) K = np.insert(K, 0, K0) for i in range(1, n): sum = 0 for j in range(0, i): sum += beta[i][j]*K[i-1] if beta[i][i] != 0: F =lambda delta: delta - f(t0+alpha[i]*h, y0+h*(delta*beta[i][i]+sum)) DF = lambda delta: np.eye(m)-Df(t0+alpha[i]*h,y0+h*(beta[i][i]*delta+sum)) sol,_ = newton(F, DF, f(t0, y0), 1e-8,1000) print(sol) else: sol = f(t0 + alpha[i]*h, y0+h*sum) K = np.insert(K, i, sol) phi = 0 for i in range(n): phi += gamma[i]*K[i] return y0 + h*phi def newton(F,DF,x0,eps,K): #---# Parameters #-------------------------------------------# #F : lambda function F(delta) #DF : our lambda function of the Jacobian of F #x0 : our initial point f(t0, y0) #eps : our tolerance for 0 #K : maximum number of steps #------------------------------------------------------------# k = 0 x = x0.copy().astype(float) delta = np.zeros([len(x0)]) Fx = F(x) while Fx.dot(Fx) > eps*eps and k<K: delta[:] = np.linalg.solve(DF(x), Fx) x[:] -= delta[:] Fx = F(x) k += 1 return [x,k] We can also use this function evolve to find the value for a given $T>t_0$ over $N$ steps so $h=(T-t_0)/N$.
def evolve(phi, f,Df, t0,y0, T,N): #---# Parameters #-------------------------------------------# #phi: our method #f(t, y): lambda function #Df(t, y): Jacobian of our function #t0: starting time #y0: initial coordinates #T: final time #N: number of steps #------------------------------------------------------------# h = T/N y = np.zeros( [N+1,len(y0)] ) y[0] = y0 t = 0 for i in range(N): y[i+1] = phi(f,Df, t,y[i], h) t = t+h return y I'm testing on the code
stepper = lambda f,Df,t0,y0,h: dirk(f,Df,t0,y0,h,aCN,bCN,gCN) y = evolve(stepper, lambda t,y: array([ -y[0], y[0]-exp(-y[2]), 1 ]), lambda t,y: array([[-1,0,0],[1,0,exp(-y[2])],[0,0,0]]), 0,array([1.,1.,0.]), 1,10) print("%d " % len(y), end="") print("%1.3e %1.3e %1.3e " % tuple(y[0]), end="") print("%1.3e %1.3e %1.3e " % tuple(y[-1]), end="") I should be getting the vector (0.3667, 1.001, 1.000) but instead I get (0.9955, 0.9955, -0.004536). I assume my code is breaking down inside the dirk function, particularly the input for Newton's Method.
