Coverage for pySDC/implementations/problem_classes/odeSystem.py: 97%

1#!/usr/bin/env python3 

2# -*- coding: utf-8 -*- 

3""" 

4Implementation of systems test problem ODEs. 

5 

6 

7Reference : 

8 

9Van der Houwen, P. J., & Sommeijer, B. P. (1991). Iterated Runge–Kutta methods 

10on parallel computers. SIAM journal on scientific and statistical computing, 

1112(5), 1000-1028. 

12""" 

13 

14import numpy as np 

15 

16from pySDC.core.errors import ProblemError 

17from pySDC.core.problem import Problem, WorkCounter 

18from pySDC.implementations.datatype_classes.mesh import mesh 

19 

20 

21class ProtheroRobinsonAutonomous(Problem): 

22 r""" 

23 Implement the Prothero-Robinson problem into autonomous form: 

24 

25 .. math:: 

26 \begin{eqnarray*} 

27 \frac{du}{dt} &=& -\frac{u^3-g(v)^3}{\epsilon} + \frac{dg}{dv}, &\quad u(0) = g(0),\\ 

28 \frac{dv}{dt} &=& 1, &\quad v(0) = 0, 

29 \end{eqnarray*} 

30 

31 with :math:`\epsilon` a stiffness parameter, that makes the problem more stiff 

32 the smaller it is (usual taken value is :math:`\epsilon=1e^{-3}`). 

33 Exact solution is given by :math:`u(t)=g(t),\;v(t)=t`, and this implementation uses 

34 :math:`g(t)=\cos(t)`. 

35 

36 Implement also the non-linear form of this problem: 

37 

38 .. math:: 

39 \frac{du}{dt} = -\frac{u^3-g(v)^3}{\epsilon} + \frac{dg}{dv}, \quad u(0) = g(0). 

40 

41 To use an other exact solution, one just have to derivate this class 

42 and overload the `g`, `dg` and `dg2` methods. For instance, 

43 to use :math:`g(t)=e^{-0.2t}`, define and use the following class: 

44 

45 >>> class MyProtheroRobinson(ProtheroRobinsonAutonomous): 

46 >>> 

47 >>> def g(self, t): 

48 >>> return np.exp(-0.2 * t) 

49 >>> 

50 >>> def dg(self, t): 

51 >>> return (-0.2) * np.exp(-0.2 * t) 

52 >>> 

53 >>> def dg2(self, t): 

54 >>> return (-0.2) ** 2 * np.exp(-0.2 * t) 

55 

56 Parameters 

57 ---------- 

58 epsilon : float, optional 

59 Stiffness parameter. The default is 1e-3. 

60 nonLinear : bool, optional 

61 Wether or not to use the non-linear form of the problem. The default is False. 

62 newton_maxiter : int, optional 

63 Maximum number of Newton iteration in solve_system. The default is 200. 

64 newton_tol : float, optional 

65 Residuum tolerance for Newton iteration in solve_system. The default is 5e-11. 

66 stop_at_nan : bool, optional 

67 Wheter to stop or not solve_system when getting NAN. The default is True. 

68 

69 References 

70 ---------- 

71 A. Prothero and A. Robinson, On the stability and accuracy of one-step methods for solving 

72 stiff systems of ordinary differential equations, Mathematics of Computation, 28 (1974), 

73 pp. 145–162. 

74 """ 

75 

76 dtype_u = mesh 

77 dtype_f = mesh 

78 

79 def __init__(self, epsilon=1e-3, nonLinear=False, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True): 

80 """Initialization routine""" 

81 nvars = 2 

82 super().__init__((nvars, None, np.dtype('float64'))) 

83 

84 self.f = self.f_NONLIN if nonLinear else self.f_LIN 

85 self.dgInv = self.dgInv_NONLIN if nonLinear else self.dgInv_LIN 

86 self._makeAttributeAndRegister( 

87 'epsilon', 'nonLinear', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True 

88 ) 

89 self.work_counters['newton'] = WorkCounter() 

90 self.work_counters['rhs'] = WorkCounter() 

91 

92 # ------------------------------------------------------------------------- 

93 # g function (analytical solution), and its first and second derivative 

94 # ------------------------------------------------------------------------- 

95 def g(self, t): 

96 r""" 

97 Function :math:`g(t) = \cos(t)`, which is the exact solution. 

98 

99 Parameters 

100 ---------- 

101 t : float 

102 Time. 

103 

104 Returns 

105 ------- 

106 float 

107 The value of :math:`g(t)`. 

108 """ 

109 return np.cos(t) 

110 

111 def dg(self, t): 

112 r""" 

113 First derivative :math:`g'(t) = -\sin(t)` of :math:`g`. 

114 

115 Parameters 

116 ---------- 

117 t : float 

118 Time. 

119 

120 Returns 

121 ------- 

122 float 

123 The value of :math:`g'(t)`. 

124 """ 

125 return -np.sin(t) 

126 

127 def dg2(self, t): 

128 r""" 

129 Second derivative :math:`g''(t) = -\cos(t)` of :math:`g`. 

130 

131 Parameters 

132 ---------- 

133 t : float 

134 Time. 

135 

136 Returns 

137 ------- 

138 float 

139 The value of :math:`g''(t)`. 

140 """ 

141 return -np.cos(t) 

142 

143 # ------------------------------------------------------------------------- 

144 # f(u,t) and Jacobian functions 

145 # ------------------------------------------------------------------------- 

146 def f(self, u, t): 

147 """ 

148 Right-hand side of the first component, replaced in ``__init__`` by ``f_LIN`` or ``f_NONLIN``. 

149 

150 Parameters 

151 ---------- 

152 u : float 

153 First solution component :math:`u`. 

154 t : float 

155 Second solution component :math:`v`, which equals the time. 

156 

157 Raises 

158 ------ 

159 NotImplementedError 

160 Always. 

161 """ 

162 raise NotImplementedError() 

163 

164 def f_LIN(self, u, t): 

165 r""" 

166 Right-hand side :math:`-(u - g(v)) / \epsilon + g'(v)` of the first component in the linear form. 

167 

168 Parameters 

169 ---------- 

170 u : float 

171 First solution component :math:`u`. 

172 t : float 

173 Second solution component :math:`v`, which equals the time. 

174 

175 Returns 

176 ------- 

177 float 

178 The right-hand side of the first component. 

179 """ 

180 return -self.epsilon ** (-1) * (u - self.g(t)) + self.dg(t) 

181 

182 def f_NONLIN(self, u, t): 

183 r""" 

184 Right-hand side :math:`-(u^3 - g(v)^3) / \epsilon + g'(v)` of the first component in the nonlinear form. 

185 

186 Parameters 

187 ---------- 

188 u : float 

189 First solution component :math:`u`. 

190 t : float 

191 Second solution component :math:`v`, which equals the time. 

192 

193 Returns 

194 ------- 

195 float 

196 The right-hand side of the first component. 

197 """ 

198 return -self.epsilon ** (-1) * (u**3 - self.g(t) ** 3) + self.dg(t) 

199 

200 def dgInv(self, u, t, dt): 

201 """ 

202 Inverse Jacobian of the Newton function, replaced in ``__init__`` by ``dgInv_LIN`` or ``dgInv_NONLIN``. 

203 

204 Parameters 

205 ---------- 

206 u : float 

207 First solution component :math:`u`. 

208 t : float 

209 Second solution component :math:`v`, which equals the time. 

210 dt : float 

211 Step size of the implicit Euler step solved by ``solve_system``. 

212 

213 Raises 

214 ------ 

215 NotImplementedError 

216 Always. 

217 """ 

218 raise NotImplementedError() 

219 

220 def dgInv_LIN(self, u, t, dt): 

221 r""" 

222 Inverse of the Jacobian of the Newton function :math:`G(u, v) = (u, v) - dt f(u, v) - rhs` of ``solve_system`` 

223 for the linear form, which does not depend on :math:`u`. 

224 

225 Parameters 

226 ---------- 

227 u : float 

228 First solution component :math:`u`. 

229 t : float 

230 Second solution component :math:`v`, which equals the time. 

231 dt : float 

232 Abbrev. for the node-to-node stepsize (or any other factor required). 

233 

234 Returns 

235 ------- 

236 np.2darray 

237 The inverse of the :math:`2 \times 2` Jacobian. 

238 """ 

239 e = self.epsilon 

240 g1, g2 = self.dg(t), self.dg2(t) 

241 return np.array([[1 / (dt / e + 1), (dt * g2 + dt * g1 / e) / (dt / e + 1)], [0, 1]]) 

242 

243 def dgInv_NONLIN(self, u, t, dt): 

244 r""" 

245 Inverse of the Jacobian of the Newton function :math:`G(u, v) = (u, v) - dt f(u, v) - rhs` of ``solve_system`` 

246 for the nonlinear form. 

247 

248 Parameters 

249 ---------- 

250 u : float 

251 First solution component :math:`u`. 

252 t : float 

253 Second solution component :math:`v`, which equals the time. 

254 dt : float 

255 Abbrev. for the node-to-node stepsize (or any other factor required). 

256 

257 Returns 

258 ------- 

259 np.2darray 

260 The inverse of the :math:`2 \times 2` Jacobian. 

261 """ 

262 e = self.epsilon 

263 g, g1, g2 = self.g(t), self.dg(t), self.dg2(t) 

264 return np.array( 

265 [[1 / (3 * dt * u**2 / e + 1), (dt * g2 + 3 * dt * g**2 * g1 / e) / (3 * dt * u**2 / e + 1)], [0, 1]] 

266 ) 

267 

268 # ------------------------------------------------------------------------- 

269 # pySDC required methods 

270 # ------------------------------------------------------------------------- 

271 def u_exact(self, t, u_init=None, t_init=None): 

272 r""" 

273 Routine to return initial conditions or exact solutions. 

274 

275 Parameters 

276 ---------- 

277 t : float 

278 Time at which the exact solution is computed. 

279 u_init : dtype_u 

280 Initial conditions for getting the exact solution. 

281 t_init : float 

282 The starting time. 

283 

284 Returns 

285 ------- 

286 u : dtype_u 

287 The exact solution. 

288 """ 

289 u = self.dtype_u(self.init) 

290 u[0] = self.g(t) 

291 u[1] = t 

292 return u 

293 

294 def eval_f(self, u, t): 

295 """ 

296 Routine to evaluate the right-hand side of the problem. 

297 

298 Parameters 

299 ---------- 

300 u : dtype_u 

301 Current values of the numerical solution. 

302 t : float 

303 Current time of the numerical solution is computed (not used here). 

304 

305 Returns 

306 ------- 

307 f : dtype_f 

308 The right-hand side of the problem (one component). 

309 """ 

310 

311 f = self.dtype_f(self.init) 

312 u, t = u 

313 f[0] = self.f(u, t) 

314 f[1] = 1 

315 self.work_counters['rhs']() 

316 return f 

317 

318 def solve_system(self, rhs, dt, u0, t): 

319 """ 

320 Simple Newton solver for the nonlinear equation 

321 

322 Parameters 

323 ---------- 

324 rhs : dtype_f 

325 Right-hand side for the nonlinear system. 

326 dt : float 

327 Abbrev. for the node-to-node stepsize (or any other factor required). 

328 u0 : dtype_u 

329 Initial guess for the iterative solver. 

330 t : float 

331 Time of the updated solution (e.g. for time-dependent BCs). 

332 

333 Returns 

334 ------- 

335 u : dtype_u 

336 The solution as mesh. 

337 """ 

338 # create new mesh object from u0 and set initial values for iteration 

339 u = self.dtype_u(u0) 

340 

341 # start newton iteration 

342 n, res = 0, np.inf 

343 while n < self.newton_maxiter: 

344 # evaluate RHS 

345 f = self.dtype_u(u) 

346 f[0] = self.f(*u) 

347 f[1] = 1 

348 

349 # form the function g with g(u) = 0 

350 g = u - dt * f - rhs 

351 

352 # if g is close to 0, then we are done 

353 res = np.linalg.norm(g, np.inf) 

354 if res < self.newton_tol or np.isnan(res): 

355 break 

356 

357 # assemble (dg/du)^{-1} 

358 dgInv = self.dgInv(u[0], u[1], dt) 

359 # newton update: u1 = u0 - g/dg 

360 u -= dgInv @ g 

361 

362 # increase iteration count and work counter 

363 n += 1 

364 self.work_counters['newton']() 

365 

366 if np.isnan(res) and self.stop_at_nan: 

367 raise ProblemError('Newton got nan after %i iterations, aborting...' % n) 

368 elif np.isnan(res): # pragma: no cover 

369 self.logger.warning('Newton got nan after %i iterations...' % n) 

370 

371 if n == self.newton_maxiter: 

372 raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res)) 

373 

374 return u 

375 

376 

377class Kaps(Problem): 

378 r""" 

379 Implement the Kaps problem: 

380 

381 .. math:: 

382 \begin{eqnarray*} 

383 \frac{du}{dt} &=& -(2+\epsilon^{-1})u + \frac{v^2}{\epsilon}, &\quad u(0) = 1,\\ 

384 \frac{dv}{dt} &=& u - v(1+v), &\quad v(0) = 1, 

385 \end{eqnarray*} 

386 

387 with :math:`\epsilon` a stiffness parameter, that makes the problem more stiff 

388 the smaller it is (usual taken value is :math:`\epsilon=1e^{-3}`). 

389 Exact solution is given by :math:`u(t)=e^{-2t},\;v(t)=e^{-t}`. 

390 

391 Parameters 

392 ---------- 

393 epsilon : float, optional 

394 Stiffness parameter. The default is 1e-3. 

395 newton_maxiter : int, optional 

396 Maximum number of Newton iteration in solve_system. The default is 200. 

397 newton_tol : float, optional 

398 Residuum tolerance for Newton iteration in solve_system. The default is 5e-11. 

399 stop_at_nan : bool, optional 

400 Wheter to stop or not solve_system when getting NAN. The default is True. 

401 

402 References 

403 ---------- 

404 Van der Houwen, P. J., & Sommeijer, B. P. (1991). Iterated Runge–Kutta methods 

405 on parallel computers. SIAM journal on scientific and statistical computing, 

406 12(5), 1000-1028. 

407 """ 

408 

409 dtype_u = mesh 

410 dtype_f = mesh 

411 

412 def __init__(self, epsilon=1e-3, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True): 

413 """Initialization routine""" 

414 nvars = 2 

415 super().__init__((nvars, None, np.dtype('float64'))) 

416 

417 self._makeAttributeAndRegister( 

418 'epsilon', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True 

419 ) 

420 self.work_counters['newton'] = WorkCounter() 

421 self.work_counters['rhs'] = WorkCounter() 

422 

423 def u_exact(self, t, u_init=None, t_init=None): 

424 r""" 

425 Routine to return initial conditions or exact solutions. 

426 

427 Parameters 

428 ---------- 

429 t : float 

430 Time at which the exact solution is computed. 

431 u_init : dtype_u 

432 Initial conditions for getting the exact solution. 

433 t_init : float 

434 The starting time. 

435 

436 Returns 

437 ------- 

438 u : dtype_u 

439 The exact solution. 

440 """ 

441 u = self.dtype_u(self.init) 

442 u[:] = [np.exp(-2 * t), np.exp(-t)] 

443 return u 

444 

445 def eval_f(self, u, t): 

446 """ 

447 Routine to evaluate the right-hand side of the problem. 

448 

449 Parameters 

450 ---------- 

451 u : dtype_u 

452 Current values of the numerical solution. 

453 t : float 

454 Current time of the numerical solution is computed (not used here). 

455 

456 Returns 

457 ------- 

458 f : dtype_f 

459 The right-hand side of the problem (one component). 

460 """ 

461 f = self.dtype_f(self.init) 

462 eps = self.epsilon 

463 x, y = u 

464 

465 f[:] = [-(2 + 1 / eps) * x + y**2 / eps, x - y * (1 + y)] 

466 self.work_counters['rhs']() 

467 return f 

468 

469 def solve_system(self, rhs, dt, u0, t): 

470 """ 

471 Simple Newton solver for the nonlinear equation 

472 

473 Parameters 

474 ---------- 

475 rhs : dtype_f 

476 Right-hand side for the nonlinear system. 

477 dt : float 

478 Abbrev. for the node-to-node stepsize (or any other factor required). 

479 u0 : dtype_u 

480 Initial guess for the iterative solver. 

481 t : float 

482 Current time (e.g. for time-dependent BCs). 

483 

484 Returns 

485 ------- 

486 u : dtype_u 

487 The solution as mesh. 

488 """ 

489 # create new mesh object from u0 and set initial values for iteration 

490 u = self.dtype_u(u0) 

491 eps = self.epsilon 

492 

493 # start newton iteration 

494 n, res = 0, np.inf 

495 while n < self.newton_maxiter: 

496 x, y = u 

497 f = np.array([-(2 + 1 / eps) * x + y**2 / eps, x - y * (1 + y)]) 

498 

499 # form the function g with g(u) = 0 

500 g = u - dt * f - rhs 

501 

502 # if g is close to 0, then we are done 

503 res = np.linalg.norm(g, np.inf) 

504 if res < self.newton_tol or np.isnan(res): 

505 break 

506 

507 # assemble (dg/du)^(-1) 

508 prefactor = 4 * dt**2 * eps * y + 2 * dt**2 * eps + dt**2 + 2 * dt * eps * y + 3 * dt * eps + dt + eps 

509 dgInv = ( 

510 1 

511 / prefactor 

512 * np.array([[2 * dt * eps * y + dt * eps + eps, 2 * dt * y], [dt * eps, 2 * dt * eps + dt + eps]]) 

513 ) 

514 

515 # newton update: u1 = u0 - g/dg 

516 u -= dgInv @ g 

517 

518 # increase iteration count and work counter 

519 n += 1 

520 self.work_counters['newton']() 

521 

522 if np.isnan(res) and self.stop_at_nan: 

523 raise ProblemError('Newton got nan after %i iterations, aborting...' % n) 

524 elif np.isnan(res): # pragma: no cover 

525 self.logger.warning('Newton got nan after %i iterations...' % n) 

526 

527 if n == self.newton_maxiter: 

528 raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res)) 

529 

530 return u 

531 

532 

533class ChemicalReaction3Var(Problem): 

534 r""" 

535 Chemical reaction with three components, modeled by the non-linear system: 

536 

537 .. math:: 

538 \frac{d{\bf u}}{dt} = 

539 \begin{pmatrix} 

540 0.013+1000u_3 & 0 & 0 \\ 

541 0 & 2500u_3 0 \\ 

542 0.013 & 0 & 1000u_1 + 2500u_2 

543 \end{pmatrix} 

544 {\bf u}, 

545 

546 with initial solution :math:`u(0)=(0.990731920827, 1.009264413846, -0.366532612659e-5)`. 

547 

548 Parameters 

549 ---------- 

550 newton_maxiter : int, optional 

551 Maximum number of Newton iteration in solve_system. The default is 200. 

552 newton_tol : float, optional 

553 Residuum tolerance for Newton iteration in solve_system. The default is 5e-11. 

554 stop_at_nan : bool, optional 

555 Wheter to stop or not solve_system when getting NAN. The default is True. 

556 

557 References 

558 ---------- 

559 Van der Houwen, P. J., & Sommeijer, B. P. (1991). Iterated Runge–Kutta methods 

560 on parallel computers. SIAM journal on scientific and statistical computing, 

561 12(5), 1000-1028. 

562 """ 

563 

564 dtype_u = mesh 

565 dtype_f = mesh 

566 

567 def __init__(self, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True): 

568 """Initialization routine""" 

569 nvars = 3 

570 u0 = (0.990731920827, 1.009264413846, -0.366532612659e-5) 

571 super().__init__((nvars, None, np.dtype('float64'))) 

572 

573 self._makeAttributeAndRegister( 

574 'u0', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True 

575 ) 

576 self.work_counters['newton'] = WorkCounter() 

577 self.work_counters['rhs'] = WorkCounter() 

578 

579 def u_exact(self, t, u_init=None, t_init=None): 

580 r""" 

581 Routine to return initial conditions or to approximate exact solution using ``SciPy``. 

582 

583 Parameters 

584 ---------- 

585 t : float 

586 Time at which the approximated exact solution is computed. 

587 u_init : pySDC.implementations.problem_classes.Lorenz.dtype_u 

588 Initial conditions for getting the exact solution. 

589 t_init : float 

590 The starting time. 

591 

592 Returns 

593 ------- 

594 me : dtype_u 

595 The approximated exact solution. 

596 """ 

597 

598 me = self.dtype_u(self.init) 

599 

600 if t > 0: 

601 

602 def eval_rhs(t, u): 

603 r""" 

604 Evaluate the right hand side, but switch the arguments for ``SciPy``. 

605 

606 Args: 

607 t (float): Time 

608 u (numpy.ndarray): Solution at time t 

609 

610 Returns: 

611 (numpy.ndarray): Right hand side evaluation 

612 """ 

613 return self.eval_f(u, t) 

614 

615 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init) 

616 else: 

617 me[:] = self.u0 

618 return me 

619 

620 def eval_f(self, u, t): 

621 """ 

622 Routine to evaluate the right-hand side of the problem. 

623 

624 Parameters 

625 ---------- 

626 u : dtype_u 

627 Current values of the numerical solution. 

628 t : float 

629 Current time of the numerical solution is computed (not used here). 

630 

631 Returns 

632 ------- 

633 f : dtype_f 

634 The right-hand side of the problem (one component). 

635 """ 

636 f = self.dtype_f(self.init) 

637 c1, c2, c3 = u 

638 

639 f[:] = -np.array([0.013 * c1 + 1000 * c3 * c1, 2500 * c3 * c2, 0.013 * c1 + 1000 * c1 * c3 + 2500 * c2 * c3]) 

640 self.work_counters['rhs']() 

641 return f 

642 

643 def solve_system(self, rhs, dt, u0, t): 

644 """ 

645 Simple Newton solver for the nonlinear equation 

646 

647 Parameters 

648 ---------- 

649 rhs : dtype_f 

650 Right-hand side for the nonlinear system. 

651 dt : float 

652 Abbrev. for the node-to-node stepsize (or any other factor required). 

653 u0 : dtype_u 

654 Initial guess for the iterative solver. 

655 t : float 

656 Current time (e.g. for time-dependent BCs). 

657 

658 Returns 

659 ------- 

660 u : dtype_u 

661 The solution as mesh. 

662 """ 

663 # create new mesh object from u0 and set initial values for iteration 

664 u = self.dtype_u(u0) 

665 

666 # start newton iteration 

667 n, res = 0, np.inf 

668 while n < self.newton_maxiter: 

669 c1, c2, c3 = u 

670 f = -np.array([0.013 * c1 + 1000 * c3 * c1, 2500 * c3 * c2, 0.013 * c1 + 1000 * c1 * c3 + 2500 * c2 * c3]) 

671 

672 # form the function g with g(u) = 0 

673 g = u - dt * f - rhs 

674 

675 # if g is close to 0, then we are done 

676 res = np.linalg.norm(g, np.inf) 

677 if res < self.newton_tol or np.isnan(res): 

678 break 

679 

680 # assemble (dg/du)^(-1) 

681 dgInv = np.array( 

682 [ 

683 [ 

684 ( 

685 2500000000.0 * c1 * c3**2 * dt**3 

686 + 32500.0 * c1 * c3 * dt**3 

687 + 3500000.0 * c1 * c3 * dt**2 

688 + 13.0 * c1 * dt**2 

689 + 1000.0 * c1 * dt 

690 + 2500000.0 * c2 * c3 * dt**2 

691 + 32.5 * c2 * dt**2 

692 + 2500.0 * c2 * dt 

693 + 2500000.0 * c3**2 * dt**2 

694 + 32.5 * c3 * dt**2 

695 + 3500.0 * c3 * dt 

696 + 0.013 * dt 

697 + 1.0 

698 ) 

699 / ( 

700 2500000000.0 * c1 * c3**2 * dt**3 

701 + 32500.0 * c1 * c3 * dt**3 

702 + 3500000.0 * c1 * c3 * dt**2 

703 + 13.0 * c1 * dt**2 

704 + 1000.0 * c1 * dt 

705 + 2500000000.0 * c2 * c3**2 * dt**3 

706 + 65000.0 * c2 * c3 * dt**3 

707 + 5000000.0 * c2 * c3 * dt**2 

708 + 0.4225 * c2 * dt**3 

709 + 65.0 * c2 * dt**2 

710 + 2500.0 * c2 * dt 

711 + 2500000000.0 * c3**3 * dt**3 

712 + 65000.0 * c3**2 * dt**3 

713 + 6000000.0 * c3**2 * dt**2 

714 + 0.4225 * c3 * dt**3 

715 + 91.0 * c3 * dt**2 

716 + 4500.0 * c3 * dt 

717 + 0.000169 * dt**2 

718 + 0.026 * dt 

719 + 1.0 

720 ), 

721 (2500000000.0 * c1 * c3**2 * dt**3 + 32500.0 * c1 * c3 * dt**3 + 2500000.0 * c1 * c3 * dt**2) 

722 / ( 

723 2500000000.0 * c1 * c3**2 * dt**3 

724 + 32500.0 * c1 * c3 * dt**3 

725 + 3500000.0 * c1 * c3 * dt**2 

726 + 13.0 * c1 * dt**2 

727 + 1000.0 * c1 * dt 

728 + 2500000000.0 * c2 * c3**2 * dt**3 

729 + 65000.0 * c2 * c3 * dt**3 

730 + 5000000.0 * c2 * c3 * dt**2 

731 + 0.4225 * c2 * dt**3 

732 + 65.0 * c2 * dt**2 

733 + 2500.0 * c2 * dt 

734 + 2500000000.0 * c3**3 * dt**3 

735 + 65000.0 * c3**2 * dt**3 

736 + 6000000.0 * c3**2 * dt**2 

737 + 0.4225 * c3 * dt**3 

738 + 91.0 * c3 * dt**2 

739 + 4500.0 * c3 * dt 

740 + 0.000169 * dt**2 

741 + 0.026 * dt 

742 + 1.0 

743 ), 

744 ( 

745 -2500000000.0 * c1 * c3**2 * dt**3 

746 - 32500.0 * c1 * c3 * dt**3 

747 - 3500000.0 * c1 * c3 * dt**2 

748 - 13.0 * c1 * dt**2 

749 - 1000.0 * c1 * dt 

750 ) 

751 / ( 

752 2500000000.0 * c1 * c3**2 * dt**3 

753 + 32500.0 * c1 * c3 * dt**3 

754 + 3500000.0 * c1 * c3 * dt**2 

755 + 13.0 * c1 * dt**2 

756 + 1000.0 * c1 * dt 

757 + 2500000000.0 * c2 * c3**2 * dt**3 

758 + 65000.0 * c2 * c3 * dt**3 

759 + 5000000.0 * c2 * c3 * dt**2 

760 + 0.4225 * c2 * dt**3 

761 + 65.0 * c2 * dt**2 

762 + 2500.0 * c2 * dt 

763 + 2500000000.0 * c3**3 * dt**3 

764 + 65000.0 * c3**2 * dt**3 

765 + 6000000.0 * c3**2 * dt**2 

766 + 0.4225 * c3 * dt**3 

767 + 91.0 * c3 * dt**2 

768 + 4500.0 * c3 * dt 

769 + 0.000169 * dt**2 

770 + 0.026 * dt 

771 + 1.0 

772 ), 

773 ], 

774 [ 

775 (6250000000.0 * c2 * c3 * dt**2 + 81250.0 * c2 * dt**2) 

776 / ( 

777 6250000000.0 * c1 * c3 * dt**2 

778 + 2500000.0 * c1 * dt 

779 + 6250000000.0 * c2 * c3 * dt**2 

780 + 81250.0 * c2 * dt**2 

781 + 6250000.0 * c2 * dt 

782 + 6250000000.0 * c3**2 * dt**2 

783 + 81250.0 * c3 * dt**2 

784 + 8750000.0 * c3 * dt 

785 + 32.5 * dt 

786 + 2500.0 

787 ), 

788 ( 

789 2500000.0 * c1 * dt 

790 + 6250000000.0 * c2 * c3 * dt**2 

791 + 81250.0 * c2 * dt**2 

792 + 6250000.0 * c2 * dt 

793 + 2500000.0 * c3 * dt 

794 + 32.5 * dt 

795 + 2500.0 

796 ) 

797 / ( 

798 6250000000.0 * c1 * c3 * dt**2 

799 + 2500000.0 * c1 * dt 

800 + 6250000000.0 * c2 * c3 * dt**2 

801 + 81250.0 * c2 * dt**2 

802 + 6250000.0 * c2 * dt 

803 + 6250000000.0 * c3**2 * dt**2 

804 + 81250.0 * c3 * dt**2 

805 + 8750000.0 * c3 * dt 

806 + 32.5 * dt 

807 + 2500.0 

808 ), 

809 (-6250000000.0 * c2 * c3 * dt**2 - 81250.0 * c2 * dt**2 - 6250000.0 * c2 * dt) 

810 / ( 

811 6250000000.0 * c1 * c3 * dt**2 

812 + 2500000.0 * c1 * dt 

813 + 6250000000.0 * c2 * c3 * dt**2 

814 + 81250.0 * c2 * dt**2 

815 + 6250000.0 * c2 * dt 

816 + 6250000000.0 * c3**2 * dt**2 

817 + 81250.0 * c3 * dt**2 

818 + 8750000.0 * c3 * dt 

819 + 32.5 * dt 

820 + 2500.0 

821 ), 

822 ], 

823 [ 

824 (-2500000.0 * c3**2 * dt**2 - 32.5 * c3 * dt**2 - 1000.0 * c3 * dt - 0.013 * dt) 

825 / ( 

826 2500000.0 * c1 * c3 * dt**2 

827 + 1000.0 * c1 * dt 

828 + 2500000.0 * c2 * c3 * dt**2 

829 + 32.5 * c2 * dt**2 

830 + 2500.0 * c2 * dt 

831 + 2500000.0 * c3**2 * dt**2 

832 + 32.5 * c3 * dt**2 

833 + 3500.0 * c3 * dt 

834 + 0.013 * dt 

835 + 1.0 

836 ), 

837 (-2500000.0 * c3**2 * dt**2 - 32.5 * c3 * dt**2 - 2500.0 * c3 * dt) 

838 / ( 

839 2500000.0 * c1 * c3 * dt**2 

840 + 1000.0 * c1 * dt 

841 + 2500000.0 * c2 * c3 * dt**2 

842 + 32.5 * c2 * dt**2 

843 + 2500.0 * c2 * dt 

844 + 2500000.0 * c3**2 * dt**2 

845 + 32.5 * c3 * dt**2 

846 + 3500.0 * c3 * dt 

847 + 0.013 * dt 

848 + 1.0 

849 ), 

850 (2500000.0 * c3**2 * dt**2 + 32.5 * c3 * dt**2 + 3500.0 * c3 * dt + 0.013 * dt + 1.0) 

851 / ( 

852 2500000.0 * c1 * c3 * dt**2 

853 + 1000.0 * c1 * dt 

854 + 2500000.0 * c2 * c3 * dt**2 

855 + 32.5 * c2 * dt**2 

856 + 2500.0 * c2 * dt 

857 + 2500000.0 * c3**2 * dt**2 

858 + 32.5 * c3 * dt**2 

859 + 3500.0 * c3 * dt 

860 + 0.013 * dt 

861 + 1.0 

862 ), 

863 ], 

864 ] 

865 ) 

866 

867 # newton update: u1 = u0 - g/dg 

868 u -= dgInv @ g 

869 

870 # increase iteration count and work counter 

871 n += 1 

872 self.work_counters['newton']() 

873 

874 if np.isnan(res) and self.stop_at_nan: 

875 raise ProblemError('Newton got nan after %i iterations, aborting...' % n) 

876 elif np.isnan(res): # pragma: no cover 

877 self.logger.warning('Newton got nan after %i iterations...' % n) 

878 

879 if n == self.newton_maxiter: 

880 raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res)) 

881 

882 return u 

883 

884 

885class JacobiElliptic(Problem): 

886 r""" 

887 Implement the Jacobi Elliptic non-linear problem: 

888 

889 .. math:: 

890 \begin{eqnarray*} 

891 \frac{du}{dt} &=& vw, &\quad u(0) = 0, \\ 

892 \frac{dv}{dt} &=& -uw, &\quad v(0) = 1, \\ 

893 \frac{dw}{dt} &=& -0.51uv, &\quad w(0) = 1. 

894 \end{eqnarray*} 

895 

896 Parameters 

897 ---------- 

898 newton_maxiter : int, optional 

899 Maximum number of Newton iteration in solve_system. The default is 200. 

900 newton_tol : float, optional