odeScalar

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.

References

A. Prothero and A. Robinson, On the stability and accuracy of one-step methods for solving stiff systems of ordinary differential equations, Mathematics of Computation, 28 (1974), pp. 145–162.

dg(t)

First derivative \(g'(t) = -\sin(t)\) of \(g\).

Parameters:

t (float) – Time.

Returns:

float – The value of \(g'(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 (dtype_f) – The right-hand side of the problem (one component).

f(u, t)

Right-hand side, replaced in __init__ by f_LIN or f_NONLIN.

Parameters:

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

• t (float) – Current time.

Raises:

NotImplementedError – Always.

f_LIN(u, t)

Right-hand side \(-(u - g(t)) / \epsilon + g'(t)\) of the linear form.

Parameters:

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

• t (float) – Current time.

Returns:

dtype_u – The right-hand side.

f_NONLIN(u, t)

Right-hand side \(-(u^3 - g(t)^3) / \epsilon + g'(t)\) of the nonlinear form.

Parameters:

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

• t (float) – Current time.

Returns:

dtype_u – The right-hand side.

g(t)

Function \(g(t) = \cos(t)\), which is the exact solution.

Parameters:

t (float) – Time.

Returns:

float – The value of \(g(t)\).

jac(u, t)

Jacobian of the right-hand side, replaced in __init__ by jac_LIN or jac_NONLIN.

Parameters:

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

• t (float) – Current time.

Raises:

NotImplementedError – Always.

jac_LIN(u, t)

Jacobian \(-1/\epsilon\) of the right-hand side of the linear form.

Parameters:

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

• t (float) – Current time.

Returns:

float – The Jacobian.

jac_NONLIN(u, t)

Jacobian \(-3 u^2 / \epsilon\) of the right-hand side of the nonlinear form.

Parameters:

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

• t (float) – Current time.

Returns:

dtype_u – The Jacobian.

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 (dtype_u) – The solution as mesh.

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 (dtype_u) – The exact solution.