Coverage for pySDC/implementations/sweeper_classes/Runge_Kutta.py: 95%

1import numpy as np 

2import logging 

3 

4from pySDC.core.Sweeper import sweeper, _Pars 

5from pySDC.core.Errors import ParameterError 

6from pySDC.core.Level import level 

7 

8 

9class ButcherTableau(object): 

10 def __init__(self, weights, nodes, matrix): 

11 """ 

12 Initialization routine to get a quadrature matrix out of a Butcher tableau 

13 

14 Args: 

15 weights (numpy.ndarray): Butcher tableau weights 

16 nodes (numpy.ndarray): Butcher tableau nodes 

17 matrix (numpy.ndarray): Butcher tableau entries 

18 """ 

19 # check if the arguments have the correct form 

20 if type(matrix) != np.ndarray: 

21 raise ParameterError('Runge-Kutta matrix needs to be supplied as a numpy array!') 

22 elif len(np.unique(matrix.shape)) != 1 or len(matrix.shape) != 2: 

23 raise ParameterError('Runge-Kutta matrix needs to be a square 2D numpy array!') 

24 

25 if type(weights) != np.ndarray: 

26 raise ParameterError('Weights need to be supplied as a numpy array!') 

27 elif len(weights.shape) != 1: 

28 raise ParameterError(f'Incompatible dimension of weights! Need 1, got {len(weights.shape)}') 

29 elif len(weights) != matrix.shape[0]: 

30 raise ParameterError(f'Incompatible number of weights! Need {matrix.shape[0]}, got {len(weights)}') 

31 

32 if type(nodes) != np.ndarray: 

33 raise ParameterError('Nodes need to be supplied as a numpy array!') 

34 elif len(nodes.shape) != 1: 

35 raise ParameterError(f'Incompatible dimension of nodes! Need 1, got {len(nodes.shape)}') 

36 elif len(nodes) != matrix.shape[0]: 

37 raise ParameterError(f'Incompatible number of nodes! Need {matrix.shape[0]}, got {len(nodes)}') 

38 

39 # Set number of nodes, left and right interval boundaries 

40 self.num_solution_stages = 1 

41 self.num_nodes = matrix.shape[0] + self.num_solution_stages 

42 self.tleft = 0.0 

43 self.tright = 1.0 

44 

45 self.nodes = np.append(np.append([0], nodes), [1]) 

46 self.weights = weights 

47 self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1]) 

48 self.Qmat[1:-1, 1:-1] = matrix 

49 self.Qmat[-1, 1:-1] = weights # this is for computing the solution to the step from the previous stages 

50 

51 self.left_is_node = True 

52 self.right_is_node = self.nodes[-1] == self.tright 

53 

54 # compute distances between the nodes 

55 if self.num_nodes > 1: 

56 self.delta_m = self.nodes[1:] - self.nodes[:-1] 

57 else: 

58 self.delta_m = np.zeros(1) 

59 self.delta_m[0] = self.nodes[0] - self.tleft 

60 

61 # check if the RK scheme is implicit 

62 self.implicit = any(matrix[i, i] != 0 for i in range(self.num_nodes - self.num_solution_stages)) 

63 

64 

65class ButcherTableauEmbedded(object): 

66 def __init__(self, weights, nodes, matrix): 

67 """ 

68 Initialization routine to get a quadrature matrix out of a Butcher tableau for embedded RK methods. 

69 

70 Be aware that the method that generates the final solution should be in the first row of the weights matrix. 

71 

72 Args: 

73 weights (numpy.ndarray): Butcher tableau weights 

74 nodes (numpy.ndarray): Butcher tableau nodes 

75 matrix (numpy.ndarray): Butcher tableau entries 

76 """ 

77 # check if the arguments have the correct form 

78 if type(matrix) != np.ndarray: 

79 raise ParameterError('Runge-Kutta matrix needs to be supplied as a numpy array!') 

80 elif len(np.unique(matrix.shape)) != 1 or len(matrix.shape) != 2: 

81 raise ParameterError('Runge-Kutta matrix needs to be a square 2D numpy array!') 

82 

83 if type(weights) != np.ndarray: 

84 raise ParameterError('Weights need to be supplied as a numpy array!') 

85 elif len(weights.shape) != 2: 

86 raise ParameterError(f'Incompatible dimension of weights! Need 2, got {len(weights.shape)}') 

87 elif len(weights[0]) != matrix.shape[0]: 

88 raise ParameterError(f'Incompatible number of weights! Need {matrix.shape[0]}, got {len(weights[0])}') 

89 

90 if type(nodes) != np.ndarray: 

91 raise ParameterError('Nodes need to be supplied as a numpy array!') 

92 elif len(nodes.shape) != 1: 

93 raise ParameterError(f'Incompatible dimension of nodes! Need 1, got {len(nodes.shape)}') 

94 elif len(nodes) != matrix.shape[0]: 

95 raise ParameterError(f'Incompatible number of nodes! Need {matrix.shape[0]}, got {len(nodes)}') 

96 

97 # Set number of nodes, left and right interval boundaries 

98 self.num_solution_stages = 2 

99 self.num_nodes = matrix.shape[0] + self.num_solution_stages 

100 self.tleft = 0.0 

101 self.tright = 1.0 

102 

103 self.nodes = np.append(np.append([0], nodes), [1, 1]) 

104 self.weights = weights 

105 self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1]) 

106 self.Qmat[1:-2, 1:-2] = matrix 

107 self.Qmat[-1, 1:-2] = weights[0] # this is for computing the higher order solution 

108 self.Qmat[-2, 1:-2] = weights[1] # this is for computing the lower order solution 

109 

110 self.left_is_node = True 

111 self.right_is_node = self.nodes[-1] == self.tright 

112 

113 # compute distances between the nodes 

114 if self.num_nodes > 1: 

115 self.delta_m = self.nodes[1:] - self.nodes[:-1] 

116 else: 

117 self.delta_m = np.zeros(1) 

118 self.delta_m[0] = self.nodes[0] - self.tleft 

119 

120 # check if the RK scheme is implicit 

121 self.implicit = any(matrix[i, i] != 0 for i in range(self.num_nodes - self.num_solution_stages)) 

122 

123 

124class RungeKutta(sweeper): 

125 nodes = None 

126 weights = None 

127 matrix = None 

128 ButcherTableauClass = ButcherTableau 

129 

130 """ 

131 Runge-Kutta scheme that fits the interface of a sweeper. 

132 Actually, the sweeper idea fits the Runge-Kutta idea when using only lower triangular rules, where solutions 

133 at the nodes are successively computed from earlier nodes. However, we only perform a single iteration of this. 

134 

135 We have two choices to realise a Runge-Kutta sweeper: We can choose Q = Q_Delta = <Butcher tableau>, but in this 

136 implementation, that would lead to a lot of wasted FLOPS from integrating with Q and then with Q_Delta and 

137 subtracting the two. For that reason, we built this new sweeper, which does not have a preconditioner. 

138 

139 This class only supports lower triangular Butcher tableaux such that the system can be solved with forward 

140 substitution. In this way, we don't get the maximum order that we could for the number of stages, but computing the 

141 stages is much cheaper. In particular, if the Butcher tableaux is strictly lower triangular, we get an explicit 

142 method, which does not require us to solve a system of equations to compute the stages. 

143 

144 Please be aware that all fundamental parameters of the Sweeper are ignored. These include 

145 

146 - num_nodes 

147 - collocation_class 

148 - initial_guess 

149 - QI 

150 

151 All of these variables are either determined by the RK rule, or are not part of an RK scheme. 

152 

153 The entries of the Butcher tableau are stored as class attributes. 

154 """ 

155 

156 def __init__(self, params): 

157 """ 

158 Initialization routine for the custom sweeper 

159 

160 Args: 

161 params: parameters for the sweeper 

162 """ 

163 # set up logger 

164 self.logger = logging.getLogger('sweeper') 

165 

166 # check if some parameters are set which only apply to actual sweepers 

167 for key in ['initial_guess', 'collocation_class', 'num_nodes']: 

168 if key in params: 

169 self.logger.warning(f'"{key}" will be ignored by Runge-Kutta sweeper') 

170 

171 # set parameters to their actual values 

172 self.coll = self.get_Butcher_tableau() 

173 params['initial_guess'] = 'zero' 

174 params['collocation_class'] = type(self.ButcherTableauClass) 

175 params['num_nodes'] = self.coll.num_nodes 

176 

177 # disable residual computation by default 

178 params['skip_residual_computation'] = params.get( 

179 'skip_residual_computation', ('IT_CHECK', 'IT_FINE', 'IT_COARSE', 'IT_UP', 'IT_DOWN') 

180 ) 

181 

182 # check if we can skip some usually unnecessary right hand side evaluations 

183 params['eval_rhs_at_right_boundary'] = params.get('eval_rhs_at_right_boundary', False) 

184 

185 self.params = _Pars(params) 

186 

187 # This will be set as soon as the sweeper is instantiated at the level 

188 self.__level = None 

189 

190 self.parallelizable = False 

191 self.QI = self.coll.Qmat 

192 

193 @classmethod 

194 def get_Q_matrix(cls): 

195 return cls.get_Butcher_tableau().Qmat 

196 

197 @classmethod 

198 def get_Butcher_tableau(cls): 

199 return cls.ButcherTableauClass(cls.weights, cls.nodes, cls.matrix) 

200 

201 @classmethod 

202 def get_update_order(cls): 

203 """ 

204 Get the order of the lower order method for doing adaptivity. Only applies to embedded methods. 

205 """ 

206 raise NotImplementedError( 

207 f"There is not an update order for RK scheme \"{cls.__name__}\" implemented. Maybe it is not an embedded scheme?" 

208 ) 

209 

210 def get_full_f(self, f): 

211 """ 

212 Get the full right hand side as a `mesh` from the right hand side 

213 

214 Args: 

215 f (dtype_f): Right hand side at a single node 

216 

217 Returns: 

218 mesh: Full right hand side as a mesh 

219 """ 

220 if type(f).__name__ in ['mesh', 'cupy_mesh']: 

221 return f 

222 elif type(f).__name__ in ['imex_mesh', 'imex_cupy_mesh']: 

223 return f.impl + f.expl 

224 elif f is None: 

225 prob = self.level.prob 

226 return self.get_full_f(prob.dtype_f(prob.init, val=0)) 

227 else: 

228 raise NotImplementedError(f'Type \"{type(f)}\" not implemented in Runge-Kutta sweeper') 

229 

230 def integrate(self): 

231 """ 

232 Integrates the right-hand side 

233 

234 Returns: 

235 list of dtype_u: containing the integral as values 

236 """ 

237 

238 # get current level and problem 

239 lvl = self.level 

240 prob = lvl.prob 

241 

242 me = [] 

243 

244 # integrate RHS over all collocation nodes 

245 for m in range(1, self.coll.num_nodes + 1): 

246 # new instance of dtype_u, initialize values with 0 

247 me.append(prob.dtype_u(prob.init, val=0.0)) 

248 for j in range(1, self.coll.num_nodes + 1): 

249 me[-1] += lvl.dt * self.coll.Qmat[m, j] * self.get_full_f(lvl.f[j]) 

250 

251 return me 

252 

253 def update_nodes(self): 

254 """ 

255 Update the u- and f-values at the collocation nodes 

256 

257 Returns: 

258 None 

259 """ 

260 

261 # get current level and problem 

262 lvl = self.level 

263 prob = lvl.prob 

264 

265 # only if the level has been touched before 

266 assert lvl.status.unlocked 

267 assert lvl.status.sweep <= 1, "RK schemes are direct solvers. Please perform only 1 iteration!" 

268 

269 # get number of collocation nodes for easier access 

270 M = self.coll.num_nodes 

271 

272 for m in range(0, M): 

273 # build rhs, consisting of the known values from above and new values from previous nodes (at k+1) 

274 rhs = prob.dtype_u(lvl.u[0]) 

275 for j in range(1, m + 1): 

276 rhs += lvl.dt * self.QI[m + 1, j] * self.get_full_f(lvl.f[j]) 

277 

278 # implicit solve with prefactor stemming from the diagonal of Qd, use previous stage as initial guess 

279 if self.coll.implicit: 

280 lvl.u[m + 1][:] = prob.solve_system( 

281 rhs, lvl.dt * self.QI[m + 1, m + 1], lvl.u[m], lvl.time + lvl.dt * self.coll.nodes[m + 1] 

282 ) 

283 else: 

284 lvl.u[m + 1][:] = rhs[:] 

285 

286 # update function values (we don't usually need to evaluate the RHS at the solution of the step) 

287 if m < M - self.coll.num_solution_stages or self.params.eval_rhs_at_right_boundary: 

288 lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], lvl.time + lvl.dt * self.coll.nodes[m + 1]) 

289 

290 # indicate presence of new values at this level 

291 lvl.status.updated = True 

292 

293 return None 

294 

295 def compute_end_point(self): 

296 """ 

297 In this Runge-Kutta implementation, the solution to the step is always stored in the last node 

298 """ 

299 self.level.uend = self.level.u[-1] 

300 

301 @property 

302 def level(self): 

303 """ 

304 Returns the current level 

305 

306 Returns: 

307 pySDC.Level.level: Current level 

308 """ 

309 return self.__level 

310 

311 @level.setter 

312 def level(self, lvl): 

313 """ 

314 Sets a reference to the current level (done in the initialization of the level) 

315 

316 Args: 

317 lvl (pySDC.Level.level): Current level 

318 """ 

319 assert isinstance(lvl, level), f"You tried to set the sweeper's level with an instance of {type(lvl)}!" 

320 if lvl.params.restol > 0: 

321 lvl.params.restol = -1 

322 self.logger.warning( 

323 'Overwriting residual tolerance with -1 because RK methods are direct and hence may not compute a residual at all!' 

324 ) 

325 

326 self.__level = lvl 

327 

328 def predict(self): 

329 """ 

330 Predictor to fill values at nodes before first sweep 

331 """ 

332 

333 # get current level and problem 

334 lvl = self.level 

335 prob = lvl.prob 

336 

337 for m in range(1, self.coll.num_nodes + 1): 

338 lvl.u[m] = prob.dtype_u(init=prob.init, val=0.0) 

339 

340 # indicate that this level is now ready for sweeps 

341 lvl.status.unlocked = True 

342 lvl.status.updated = True 

343 

344 

345class RungeKuttaIMEX(RungeKutta): 

346 """ 

347 Implicit-explicit split Runge Kutta base class. Only supports methods that share the nodes and weights. 

348 """ 

349 

350 matrix_explicit = None 

351 ButcherTableauClass_explicit = ButcherTableau 

352 

353 def __init__(self, params): 

354 """ 

355 Initialization routine 

356 

357 Args: 

358 params: parameters for the sweeper 

359 """ 

360 super().__init__(params) 

361 self.coll_explicit = self.get_Butcher_tableau_explicit() 

362 self.QE = self.coll_explicit.Qmat 

363 

364 def predict(self): 

365 """ 

366 Predictor to fill values at nodes before first sweep 

367 """ 

368 

369 # get current level and problem 

370 lvl = self.level 

371 prob = lvl.prob 

372 

373 for m in range(1, self.coll.num_nodes + 1): 

374 lvl.u[m] = prob.dtype_u(init=prob.init, val=0.0) 

375 lvl.f[m] = prob.dtype_f(init=prob.init, val=0.0) 

376 

377 # indicate that this level is now ready for sweeps 

378 lvl.status.unlocked = True 

379 lvl.status.updated = True 

380 

381 @classmethod 

382 def get_Butcher_tableau_explicit(cls): 

383 return cls.ButcherTableauClass_explicit(cls.weights, cls.nodes, cls.matrix_explicit) 

384 

385 def integrate(self): 

386 """ 

387 Integrates the right-hand side 

388 

389 Returns: 

390 list of dtype_u: containing the integral as values 

391 """ 

392 

393 # get current level and problem 

394 lvl = self.level 

395 prob = lvl.prob 

396 

397 me = [] 

398 

399 # integrate RHS over all collocation nodes 

400 for m in range(1, self.coll.num_nodes + 1): 

401 # new instance of dtype_u, initialize values with 0 

402 me.append(prob.dtype_u(prob.init, val=0.0)) 

403 for j in range(1, self.coll.num_nodes + 1): 

404 me[-1] += lvl.dt * ( 

405 self.coll.Qmat[m, j] * lvl.f[j].impl + self.coll_explicit.Qmat[m, j] * lvl.f[j].expl 

406 ) 

407 

408 return me 

409 

410 def update_nodes(self): 

411 """ 

412 Update the u- and f-values at the collocation nodes 

413 

414 Returns: 

415 None 

416 """ 

417 

418 # get current level and problem 

419 lvl = self.level 

420 prob = lvl.prob 

421 

422 # only if the level has been touched before 

423 assert lvl.status.unlocked 

424 assert lvl.status.sweep <= 1, "RK schemes are direct solvers. Please perform only 1 iteration!" 

425 

426 # get number of collocation nodes for easier access 

427 M = self.coll.num_nodes 

428 

429 for m in range(0, M): 

430 # build rhs, consisting of the known values from above and new values from previous nodes (at k+1) 

431 rhs = lvl.u[0] 

432 for j in range(1, m + 1): 

433 rhs += lvl.dt * (self.QI[m + 1, j] * lvl.f[j].impl + self.QE[m + 1, j] * lvl.f[j].expl) 

434 

435 # implicit solve with prefactor stemming from the diagonal of Qd, use previous stage as initial guess 

436 lvl.u[m + 1][:] = prob.solve_system( 

437 rhs, lvl.dt * self.QI[m + 1, m + 1], lvl.u[m], lvl.time + lvl.dt * self.coll.nodes[m + 1] 

438 ) 

439 

440 # update function values (we don't usually need to evaluate the RHS at the solution of the step) 

441 if m < M - self.coll.num_solution_stages or self.params.eval_rhs_at_right_boundary: 

442 lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], lvl.time + lvl.dt * self.coll.nodes[m + 1]) 

443 

444 # indicate presence of new values at this level 

445 lvl.status.updated = True 

446 

447 return None 

448 

449 

450class ForwardEuler(RungeKutta): 

451 """ 

452 Forward Euler. Still a classic. 

453 

454 Not very stable first order method. 

455 """ 

456 

457 nodes = np.array([0.0]) 

458 weights = np.array([1.0]) 

459 matrix = np.array( 

460 [ 

461 [0.0], 

462 ] 

463 ) 

464 

465 

466class BackwardEuler(RungeKutta): 

467 """ 

468 Backward Euler. A favorite among true connoisseurs of the heat equation. 

469 

470 A-stable first order method. 

471 """ 

472 

473 nodes = np.array([1.0]) 

474 weights = np.array([1.0]) 

475 matrix = np.array( 

476 [ 

477 [1.0], 

478 ] 

479 ) 

480 

481 

482class CrankNicholson(RungeKutta): 

483 """ 

484 Implicit Runge-Kutta method of second order, A-stable. 

485 """ 

486 

487 nodes = np.array([0, 1]) 

488 weights = np.array([0.5, 0.5]) 

489 matrix = np.zeros((2, 2)) 

490 matrix[1, 0] = 0.5 

491 matrix[1, 1] = 0.5 

492 

493 

494class ExplicitMidpointMethod(RungeKutta): 

495 """ 

496 Explicit Runge-Kutta method of second order. 

497 """ 

498 

499 nodes = np.array([0, 0.5]) 

500 weights = np.array([0, 1]) 

501 matrix = np.zeros((2, 2)) 

502 matrix[1, 0] = 0.5 

503 

504 

505class ImplicitMidpointMethod(RungeKutta): 

506 """ 

507 Implicit Runge-Kutta method of second order. 

508 """ 

509 

510 nodes = np.array([0.5]) 

511 weights = np.array([1]) 

512 matrix = np.zeros((1, 1)) 

513 matrix[0, 0] = 1.0 / 2.0 

514 

515 

516class RK4(RungeKutta): 

517 """ 

518 Explicit Runge-Kutta of fourth order: Everybody's darling. 

519 """ 

520 

521 nodes = np.array([0, 0.5, 0.5, 1]) 

522 weights = np.array([1.0, 2.0, 2.0, 1.0]) / 6.0 

523 matrix = np.zeros((4, 4)) 

524 matrix[1, 0] = 0.5 

525 matrix[2, 1] = 0.5 

526 matrix[3, 2] = 1.0 

527 

528 

529class Heun_Euler(RungeKutta): 

530 """ 

531 Second order explicit embedded Runge-Kutta method. 

532 """ 

533 

534 nodes = np.array([0, 1]) 

535 weights = np.array([[0.5, 0.5], [1, 0]]) 

536 matrix = np.zeros((2, 2)) 

537 matrix[1, 0] = 1 

538 ButcherTableauClass = ButcherTableauEmbedded 

539 

540 @classmethod 

541 def get_update_order(cls): 

542 return 2 

543 

544 

545class Cash_Karp(RungeKutta): 

546 """ 

547 Fifth order explicit embedded Runge-Kutta. See [here](https://doi.org/10.1145/79505.79507). 

548 """ 

549 

550 nodes = np.array([0, 0.2, 0.3, 0.6, 1.0, 7.0 / 8.0]) 

551 weights = np.array( 

552 [ 

553 [37.0 / 378.0, 0.0, 250.0 / 621.0, 125.0 / 594.0, 0.0, 512.0 / 1771.0], 

554 [2825.0 / 27648.0, 0.0, 18575.0 / 48384.0, 13525.0 / 55296.0, 277.0 / 14336.0, 1.0 / 4.0], 

555 ] 

556 ) 

557 matrix = np.zeros((6, 6)) 

558 matrix[1, 0] = 1.0 / 5.0 

559 matrix[2, :2] = [3.0 / 40.0, 9.0 / 40.0] 

560 matrix[3, :3] = [0.3, -0.9, 1.2] 

561 matrix[4, :4] = [-11.0 / 54.0, 5.0 / 2.0, -70.0 / 27.0, 35.0 / 27.0] 

562 matrix[5, :5] = [1631.0 / 55296.0, 175.0 / 512.0, 575.0 / 13824.0, 44275.0 / 110592.0, 253.0 / 4096.0] 

563 ButcherTableauClass = ButcherTableauEmbedded 

564 

565 @classmethod 

566 def get_update_order(cls): 

567 return 5 

568 

569 

570class DIRK43(RungeKutta): 

571 """ 

572 Embedded A-stable diagonally implicit RK pair of order 3 and 4. 

573 

574 Taken from [here](https://doi.org/10.1007/BF01934920). 

575 """ 

576 

577 nodes = np.array([5.0 / 6.0, 10.0 / 39.0, 0, 1.0 / 6.0]) 

578 weights = np.array( 

579 [[61.0 / 150.0, 2197.0 / 2100.0, 19.0 / 100.0, -9.0 / 14.0], [32.0 / 75.0, 169.0 / 300.0, 1.0 / 100.0, 0.0]] 

580 ) 

581 matrix = np.zeros((4, 4)) 

582 matrix[0, 0] = 5.0 / 6.0 

583 matrix[1, :2] = [-15.0 / 26.0, 5.0 / 6.0] 

584 matrix[2, :3] = [215.0 / 54.0, -130.0 / 27.0, 5.0 / 6.0] 

585 matrix[3, :] = [4007.0 / 6075.0, -31031.0 / 24300.0, -133.0 / 2700.0, 5.0 / 6.0] 

586 ButcherTableauClass = ButcherTableauEmbedded 

587 

588 @classmethod 

589 def get_update_order(cls): 

590 return 4 

591 

592 

593class ESDIRK53(RungeKutta): 

594 """ 

595 A-stable embedded RK pair of orders 5 and 3, ESDIRK5(3)6L[2]SA. 

596 Taken from [here](https://ntrs.nasa.gov/citations/20160005923) 

597 """ 

598 

599 nodes = np.array( 

600 [0, 4024571134387.0 / 7237035672548.0, 14228244952610.0 / 13832614967709.0, 1.0 / 10.0, 3.0 / 50.0, 1.0] 

601 ) 

602 matrix = np.zeros((6, 6)) 

603 matrix[1, :2] = [3282482714977.0 / 11805205429139.0, 3282482714977.0 / 11805205429139.0] 

604 matrix[2, :3] = [ 

605 606638434273.0 / 1934588254988, 

606 2719561380667.0 / 6223645057524, 

607 3282482714977.0 / 11805205429139.0, 

608 ] 

609 matrix[3, :4] = [ 

610 -651839358321.0 / 6893317340882, 

611 -1510159624805.0 / 11312503783159, 

612 235043282255.0 / 4700683032009.0, 

613 3282482714977.0 / 11805205429139.0, 

614 ] 

615 matrix[4, :5] = [ 

616 -5266892529762.0 / 23715740857879, 

617 -1007523679375.0 / 10375683364751, 

618 521543607658.0 / 16698046240053.0, 

619 514935039541.0 / 7366641897523.0, 

620 3282482714977.0 / 11805205429139.0, 

621 ] 

622 matrix[5, :] = [ 

623 -6225479754948.0 / 6925873918471, 

624 6894665360202.0 / 11185215031699, 

625 -2508324082331.0 / 20512393166649, 

626 -7289596211309.0 / 4653106810017.0, 

627 39811658682819.0 / 14781729060964.0, 

628 3282482714977.0 / 11805205429139, 

629 ] 

630 

631 weights = np.array( 

632 [ 

633 [ 

634 -6225479754948.0 / 6925873918471, 

635 6894665360202.0 / 11185215031699.0, 

636 -2508324082331.0 / 20512393166649, 

637 -7289596211309.0 / 4653106810017, 

638 39811658682819.0 / 14781729060964.0, 

639 3282482714977.0 / 11805205429139, 

640 ], 

641 [ 

642 -2512930284403.0 / 5616797563683, 

643 5849584892053.0 / 8244045029872, 

644 -718651703996.0 / 6000050726475.0, 

645 -18982822128277.0 / 13735826808854.0, 

646 23127941173280.0 / 11608435116569.0, 

647 2847520232427.0 / 11515777524847.0, 

648 ], 

649 ] 

650 ) 

651 ButcherTableauClass = ButcherTableauEmbedded 

652 

653 @classmethod 

654 def get_update_order(cls): 

655 return 4 

656 

657 

658class ESDIRK43(RungeKutta): 

659 """ 

660 A-stable embedded RK pair of orders 4 and 3, ESDIRK4(3)6L[2]SA. 

661 Taken from [here](https://ntrs.nasa.gov/citations/20160005923) 

662 """ 

663 

664 s2 = 2**0.5 

665 

666 nodes = np.array([0, 1 / 2, (2 - 2**0.5) / 4, 5 / 8, 26 / 25, 1.0]) 

667 matrix = np.zeros((6, 6)) 

668 matrix[1, :2] = [1 / 4, 1 / 4] 

669 matrix[2, :3] = [ 

670 (1 - 2**0.5) / 8, 

671 (1 - 2**0.5) / 8, 

672 1 / 4, 

673 ] 

674 matrix[3, :4] = [ 

675 (5 - 7 * s2) / 64, 

676 (5 - 7 * s2) / 64, 

677 7 * (1 + s2) / 32, 

678 1 / 4, 

679 ] 

680 matrix[4, :5] = [ 

681 (-13796 - 54539 * s2) / 125000, 

682 (-13796 - 54539 * s2) / 125000, 

683 (506605 + 132109 * s2) / 437500, 

684 166 * (-97 + 376 * s2) / 109375, 

685 1 / 4, 

686 ] 

687 matrix[5, :] = [ 

688 (1181 - 987 * s2) / 13782, 

689 (1181 - 987 * s2) / 13782, 

690 47 * (-267 + 1783 * s2) / 273343, 

691 -16 * (-22922 + 3525 * s2) / 571953, 

692 -15625 * (97 + 376 * s2) / 90749876, 

693 1 / 4, 

694 ] 

695 

696 weights = np.array( 

697 [ 

698 [ 

699 (1181 - 987 * s2) / 13782, 

700 (1181 - 987 * s2) / 13782, 

701 47 * (-267 + 1783 * s2) / 273343, 

702 -16 * (-22922 + 3525 * s2) / 571953, 

703 -15625 * (97 + 376 * s2) / 90749876, 

704 1 / 4, 

705 ], 

706 [ 

707 -480923228411.0 / 4982971448372, 

708 -480923228411.0 / 4982971448372, 

709 6709447293961.0 / 12833189095359, 

710 3513175791894.0 / 6748737351361.0, 

711 -498863281070.0 / 6042575550617.0, 

712 2077005547802.0 / 8945017530137.0, 

713 ], 

714 ] 

715 ) 

716 ButcherTableauClass = ButcherTableauEmbedded 

717 

718 @classmethod 

719 def get_update_order(cls): 

720 return 4 

721 

722 

723class ARK548L2SAERK(RungeKutta): 

724 """ 

725 Explicit part of the ARK54 scheme. 

726 """ 

727 

728 ButcherTableauClass = ButcherTableauEmbedded 

729 weights = np.array( 

730 [ 

731 [ 

732 -872700587467.0 / 9133579230613.0, 

733 0.0, 

734 0.0, 

735 22348218063261.0 / 9555858737531.0, 

736 -1143369518992.0 / 8141816002931.0, 

737 -39379526789629.0 / 19018526304540.0, 

738 32727382324388.0 / 42900044865799.0, 

739 41.0 / 200.0, 

740 ], 

741 [ 

742 -975461918565.0 / 9796059967033.0, 

743 0.0, 

744 0.0, 

745 78070527104295.0 / 32432590147079.0, 

746 -548382580838.0 / 3424219808633.0, 

747 -33438840321285.0 / 15594753105479.0, 

748 3629800801594.0 / 4656183773603.0, 

749 4035322873751.0 / 18575991585200.0, 

750 ], 

751 ] 

752 ) 

753 

754 nodes = np.array( 

755 [ 

756 0, 

757 41.0 / 100.0, 

758 2935347310677.0 / 11292855782101.0, 

759 1426016391358.0 / 7196633302097.0, 

760 92.0 / 100.0, 

761 24.0 / 100.0, 

762 3.0 / 5.0, 

763 1.0, 

764 ] 

765 ) 

766 

767 matrix = np.zeros((8, 8)) 

768 matrix[1, 0] = 41.0 / 100.0 

769 matrix[2, :2] = [367902744464.0 / 2072280473677.0, 677623207551.0 / 8224143866563.0] 

770 matrix[3, :3] = [1268023523408.0 / 10340822734521.0, 0.0, 1029933939417.0 / 13636558850479.0] 

771 matrix[4, :4] = [ 

772 14463281900351.0 / 6315353703477.0, 

773 0.0, 

774 66114435211212.0 / 5879490589093.0, 

775 -54053170152839.0 / 4284798021562.0, 

776 ] 

777 matrix[5, :5] = [ 

778 14090043504691.0 / 34967701212078.0, 

779 0.0, 

780 15191511035443.0 / 11219624916014.0,