Source code for naginterfaces.library.examples.ode.ivp_bdf_zero_simple_ex

#!/usr/bin/env python3
"``naginterfaces.library.ode.ivp_bdf_zero_simple`` Python Example."

# NAG Copyright 2018-2019.

# pylint: disable=invalid-name,too-many-locals

from naginterfaces.library import ode

def main():
    """
    Example for :func:`naginterfaces.library.ode.ivp_bdf_zero_simple`.

    Solves a stiff system of ordinary differential equations
    by the Backward Differentiation Formulae.

    Demonstrates termination by root-finding condition and intermediate
    solution output.

    >>> main()
    naginterfaces.library.ode.ivp_bdf_zero_simple Python Example Results.
    Solve stiff Robertson problem until a stopping criterion is met.
          x                       y
        0.00      1.00000      0.00000      0.00000
        2.00      0.94161      0.00003      0.05837
        4.00      0.90551      0.00002      0.09446
    Stopping criterion met at x =   4.377
    Solution is       0.90000      0.00002      0.09998
    """

    print(
        'naginterfaces.library.ode.ivp_bdf_zero_simple Python Example Results.'
    )
    print('Solve stiff Robertson problem until a stopping criterion is met.')

    # The ODE system definition:
    cb_fcn = lambda x, y: [
        -alpha*y[0] + beta*y[1]*y[2],
        alpha*y[0] - beta*y[1]*y[2] - gamma*y[1]*y[1],
        gamma*y[1]*y[1],
    ]

    # The Jacobian of the ODE system:
    cb_pederv = lambda x, y: [
        [-alpha, beta*y[2], beta*y[1]],
        [alpha, -beta*y[2] - 2.0*gamma*y[1], -beta*y[1]],
        [0., 2.*gamma*y[1], 0.],
    ]

    # Stopping criterion: y[0] attains value y0_end:
    cb_g = lambda x, y: y[0] - y0_end

    def cb_output(xsol, y):
        """Intermediate output: at steps of h."""
        print(
            '{:8.2f}'.format(xsol) +
            ''.join(['{:13.5f}']*len(y)).format(*y)
        )
        xsol += h
        # Last output point is x_end:
        if abs(xsol-x_end) < 0.25*h:
            xsol = x_end
        return xsol

    # Problem parameters:
    alpha = 0.04
    beta = 10000.
    gamma = 30000000.

    # Integration limits, stopping value and output steps:
    x_init = 0.
    x_end = 10.
    y0_end = 0.9
    h = 2.

    # Initial profile:
    y_init = [1., 0., 0.]

    # Algorithmic parameters:
    tol = 0.0001
    relabs = 'Default'

    # Output header:
    print('      x                       y')

    # Solve ODE system and print solution at intermediate points:
    ode_soln = ode.ivp_bdf_zero_simple(
        x_init, x_end, y_init, cb_fcn, tol, relabs,
        pederv=cb_pederv, output=cb_output, g=cb_g,
    )

    print('Stopping criterion met at x = {:7.3f}'.format(ode_soln.x))
    print(
        'Solution is ' +
        ''.join(['{:13.5f}']*len(ode_soln.y)).format(*ode_soln.y)
    )

if __name__ == '__main__':
    import doctest
    import sys
    sys.exit(
        doctest.testmod(
            None, verbose=True, report=False,
            optionflags=doctest.REPORT_NDIFF,
        ).failed
    )