#!/usr/bin/env python3
"``naginterfaces.library.ode.ivp_stiff_exp_fulljac`` Python Example."
# NAG Copyright 2017-2019.
# pylint: disable=invalid-name,too-many-locals
from naginterfaces.base import utils
from naginterfaces.library import ode
[docs]def main():
"""
Example for :func:`naginterfaces.library.ode.ivp_stiff_exp_fulljac`.
Explicit ODEs, stiff initial value problem, full Jacobian.
>>> main()
naginterfaces.library.ode.ivp_stiff_exp_fulljac Python Example Results.
Integrate the stiff Robertson problem.
naginterfaces.base.ode.ivp_stiff_exp_fulljac: INITIAL STEP SIZE IS ...
...
At time 10.0 the corresponding y is (0.8414, 0.0000, 0.1586).
"""
print(
'naginterfaces.library.ode.ivp_stiff_exp_fulljac '
'Python Example Results.'
)
print('Integrate the stiff Robertson problem.')
# The number of ODEs to solve:
neq = 3
neqmax = 3
# Parameters for the BDF setup routine:
maxord = 5
sdysav = maxord + 1
method = 'NewtonMod'
petzld = False
con = [0.]*6
tcrit = 10.
hmin = 1.e-10
hmax = 10.
h0 = 0.
maxstp = 200
mxhnil = 5
norm = 'AveragedL2'
con, comm = ode.ivp_stiff_bdf(
neqmax, sdysav, maxord, method, petzld, con, tcrit, hmin, hmax, h0,
maxstp, mxhnil, norm,
)
# The linear algebra setup:
jceval = 'Numerical'
nwkjac = neqmax*(neqmax+1)
ode.ivp_stiff_fulljac_setup(neq, neqmax, jceval, nwkjac, comm)
# The initial values:
y = [1., 0., 0.]
# The initial and next integration points:
t = 0.
tout = tcrit
# Tolerance settings:
itol = 1
rtol = [1.e-4]
atol = [1.e-7]
# Parameters for the Robertson problem:
alpha = 0.04
beta = 1.e4
gamma = 3.e7
def cb_fcn(t, y, ires):
"""The function for the derivative vector."""
return (
[
-alpha*y[0] + beta*y[1]*y[2],
alpha*y[0] - beta*y[1]*y[2] - gamma*y[1]**2,
gamma*y[1]**2,
],
ires
)
# Normal computation without overshooting:
itask = 4
# Enable some iteration tracing, to the NAG Engine
# advisory unit (associated with sys.stdout by default):
itrace = 1
iom = utils.FileObjManager()
# Jacobian computed numerically, default jac (None) is
# acceptable:
ode_soln = ode.ivp_stiff_exp_fulljac(
t, tout, y, comm, rtol, atol, itol, cb_fcn, itask, itrace,
io_manager=iom,
)
print(
(
'At time {:.1f} the corresponding y is (' +
', '.join(['{:.4f}']*len(ode_soln.y)) +
').'
).format(ode_soln.t, *ode_soln.y)
)
if __name__ == '__main__':
import doctest
import sys
sys.exit(
doctest.testmod(
None, verbose=True, report=False,
optionflags=doctest.ELLIPSIS | doctest.REPORT_NDIFF,
).failed
)
