implementations.problem_classes.odeScalar module

Implementation of scalar test problem ODEs.

Reference :

Van der Houwen, P. J., & Sommeijer, B. P. (1991). Iterated Runge–Kutta methods on parallel computers. SIAM journal on scientific and statistical computing, 12(5), 1000-1028.

class ProtheroRobinson(epsilon=0.001, nonLinear=False, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True)

Bases: Problem

Implement the Prothero-Robinson problem:

\[\frac{du}{dt} = -\frac{u-g(t)}{\epsilon} + \frac{dg}{dt}, \quad u(0) = g(0).,\]

with \(\epsilon\) a stiffness parameter, that makes the problem more stiff the smaller it is (usual taken value is \(\epsilon=1e^{-3}\)). Exact solution is given by \(u(t)=g(t)\), and this implementation uses \(g(t)=\cos(t)\).

Implement also the non-linear form of this problem:

\[\frac{du}{dt} = -\frac{u^3-g(t)^3}{\epsilon} + \frac{dg}{dt}, \quad u(0) = g(0).\]

To use an other exact solution, one just have to derivate this class and overload the g and dg methods. For instance, to use \(g(t)=e^{-0.2*t}\), define and use the following class:

>>> class MyProtheroRobinson(ProtheroRobinson):
>>>
>>>     def g(self, t):
>>>         return np.exp(-0.2 * t)
>>>
>>>     def dg(self, t):
>>>         return (-0.2) * np.exp(-0.2 * t)

Parameters:

• epsilon (float, optional) – Stiffness parameter. The default is 1e-3.

• nonLinear (bool, optional) – Wether or not to use the non-linear form of the problem. The default is False.

• newton_maxiter (int, optional) – Maximum number of Newton iteration in solve_system. The default is 200.

• newton_tol (float, optional) – Residuum tolerance for Newton iteration in solve_system. The default is 5e-11.

• stop_at_nan (bool, optional) – Wheter to stop or not solve_system when getting NAN. The default is True.

• Reference

• ---------

• Robinson (A. Prothero and A.)

• solving (On the stability and accuracy of one-step methods for)

• equations (stiff systems of ordinary differential)

• Computation (Mathematics of)

• (1974) (28)

:param : :param pp. 145–162.:

dg(t)

dtype_f

alias of mesh

dtype_u

alias of mesh

eval_f(u, t)

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

Parameters:

• u (dtype_u) – Current values of the numerical solution.

• t (float) – Current time of the numerical solution is computed (not used here).

Returns:

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

Return type:

dtype_f

f(u, t)

f_LIN(u, t)

f_NONLIN(u, t)

g(t)

jac(u, t)

jac_LIN(u, t)

jac_NONLIN(u, t)

solve_system(rhs, dt, u0, t)

Simple Newton solver for the nonlinear equation

Parameters:

• rhs (dtype_f) – Right-hand side for the nonlinear system.

• dt (float) – Abbrev. for the node-to-node stepsize (or any other factor required).

• u0 (dtype_u) – Initial guess for the iterative solver.

• t (float) – Time of the updated solution (e.g. for time-dependent BCs).

Returns:

u – The solution as mesh.

Return type:

dtype_u

u_exact(t, u_init=None, t_init=None)

Routine to return initial conditions or exact solution.

Parameters:

• t (float) – Time at which the exact solution is computed.

• u_init (dtype_u) – Initial conditions for getting the exact solution.

• t_init (float) – The starting time.

Returns:

u – The exact solution.

Return type:

dtype_u