Coverage for pySDC/implementations/problem_classes/Van_der_Pol_implicit.py: 98%
64 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-09-29 12:50 +0000
1import numpy as np
3from pySDC.core.errors import ProblemError
4from pySDC.core.problem import Problem, WorkCounter
5from pySDC.implementations.datatype_classes.mesh import mesh
8# noinspection PyUnusedLocal
9class vanderpol(Problem):
10 r"""
11 Stiff Van der Pol oscillator as a system of two first-order ODEs, fully implicit with Newton.
13 It is given by the equation
15 .. math::
16 \frac{d^2 u(t)}{d t^2} - \mu (1 - u(t)^2) \frac{d u(t)}{dt} + u(t) = 0.
18 Parameters
19 ----------
20 u0 : sequence of array_like, optional
21 Initial condition.
22 mu : float, optional
23 Stiff parameter :math:`\mu`.
24 newton_maxiter : int, optional
25 Maximum number of iterations for Newton's method to terminate.
26 newton_tol : float, optional
27 Tolerance for Newton to terminate.
28 stop_at_nan : bool, optional
29 Indicate whether Newton's method should stop if ``nan`` values arise.
30 crash_at_maxiter : bool, optional
31 Indicates whether Newton's method should stop if maximum number of iterations
32 ``newton_maxiter`` is reached.
33 relative_tolerance : bool, optional
34 Use a relative or absolute tolerance for the Newton solver
36 Attributes
37 ----------
38 work_counters : WorkCounter
39 Counts different things, here: Number of evaluations of the right-hand side in ``eval_f``
40 and number of Newton calls in each Newton iterations are counted.
41 """
43 dtype_u = mesh
44 dtype_f = mesh
46 def __init__(
47 self,
48 u0=None,
49 mu=5.0,
50 newton_maxiter=100,
51 newton_tol=1e-9,
52 stop_at_nan=True,
53 crash_at_maxiter=True,
54 relative_tolerance=False,
55 ):
56 """Initialization routine"""
57 nvars = 2
59 if u0 is None:
60 u0 = [2.0, 0.0]
62 super().__init__((nvars, None, np.dtype('float64')))
63 self._makeAttributeAndRegister('nvars', 'u0', localVars=locals(), readOnly=True)
64 self._makeAttributeAndRegister(
65 'mu',
66 'newton_maxiter',
67 'newton_tol',
68 'stop_at_nan',
69 'crash_at_maxiter',
70 'relative_tolerance',
71 localVars=locals(),
72 )
73 self.work_counters['newton'] = WorkCounter()
74 self.work_counters['jacobian_solves'] = WorkCounter()
75 self.work_counters['rhs'] = WorkCounter()
77 def u_exact(self, t, u_init=None, t_init=None):
78 r"""
79 Routine to approximate the exact solution at time t by ``SciPy`` or give initial conditions when called at :math:`t=0`.
81 Parameters
82 ----------
83 t : float
84 Current time.
85 u_init : pySDC.problem.vanderpol.dtype_u
86 Initial conditions for getting the exact solution.
87 t_init : float
88 The starting time.
90 Returns
91 -------
92 me : dtype_u
93 Approximate exact solution.
94 """
96 me = self.dtype_u(self.init)
98 if t > 0.0:
100 def eval_rhs(t, u):
101 return self.eval_f(u, t)
103 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
104 else:
105 me[:] = self.u0
106 return me
108 def eval_f(self, u, t):
109 """
110 Routine to compute the right-hand side for both components simultaneously.
112 Parameters
113 ----------
114 u : dtype_u
115 Current values of the numerical solution.
116 t : float
117 Current time at which the numerical solution is computed (not used here).
119 Returns
120 -------
121 f : dtype_f
122 The right-hand side (contains 2 components).
123 """
125 x1 = u[0]
126 x2 = u[1]
127 f = self.f_init
128 f[0] = x2
129 f[1] = self.mu * (1 - x1**2) * x2 - x1
130 self.work_counters['rhs']()
131 return f
133 def solve_system(self, rhs, dt, u0, t):
134 """
135 Simple Newton solver for the nonlinear system.
137 Parameters
138 ----------
139 rhs : dtype_f
140 Right-hand side for the nonlinear system.
141 dt : float
142 Abbrev. for the node-to-node stepsize (or any other factor required).
143 u0 : dtype_u
144 Initial guess for the iterative solver.
145 t : float
146 Current time (e.g. for time-dependent BCs).
148 Returns
149 -------
150 u : dtype_u
151 The solution u.
152 """
154 mu = self.mu
156 # create new mesh object from u0 and set initial values for iteration
157 u = self.dtype_u(u0)
158 x1 = u[0]
159 x2 = u[1]
161 # start newton iteration
162 n = 0
163 res = 99
164 while n < self.newton_maxiter:
165 # form the function g with g(u) = 0
166 g = np.array([x1 - dt * x2 - rhs[0], x2 - dt * (mu * (1 - x1**2) * x2 - x1) - rhs[1]])
168 # if g is close to 0, then we are done
169 res = np.linalg.norm(g, np.inf) / (abs(u) if self.relative_tolerance else 1.0)
170 if res < self.newton_tol or np.isnan(res):
171 break
173 u -= self.solve_jacobian(g, dt, u)
175 # set new values and increase iteration count
176 x1 = u[0]
177 x2 = u[1]
178 n += 1
179 self.work_counters['newton']()
181 if np.isnan(res) and self.stop_at_nan:
182 self.logger.warning('Newton got nan after %i iterations...' % n)
183 raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
184 elif np.isnan(res):
185 self.logger.warning('Newton got nan after %i iterations...' % n)
187 if n == self.newton_maxiter and self.crash_at_maxiter:
188 raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
190 return u
192 def solve_jacobian(self, rhs, dt, u, **kwargs):
193 r"""
194 Solve the linear system with the Jacobian of the Newton function :math:`g(u) = u - dt f(u) - rhs` of
195 ``solve_system``, by applying the analytically computed inverse of the :math:`2 \times 2` Jacobian at ``u``.
197 Parameters
198 ----------
199 rhs : np.1darray
200 Right-hand side of the linear system, i.e., the Newton residual.
201 dt : float
202 Abbrev. for the node-to-node stepsize (or any other factor required).
203 u : dtype_u
204 Current Newton iterate, at which the Jacobian is evaluated.
205 **kwargs
206 Not used.
208 Returns
209 -------
210 du : np.1darray
211 The Newton update.
212 """
213 mu = self.mu
214 u1 = u[0]
215 u2 = u[1]
217 # assemble prefactor
218 c = 1.0 / (-2 * dt**2 * mu * u1 * u2 - dt**2 - 1 + dt * mu * (1 - u1**2))
219 # assemble (dg/du)^-1
220 dg = c * np.array([[dt * mu * (1 - u1**2) - 1, -dt], [2 * dt * mu * u1 * u2 + dt, -1]])
222 self.work_counters['jacobian_solves']()
223 return np.dot(dg, rhs)