In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import sympy
from pyodesys.symbolic import SymbolicSys
from pyodesys.tests._robertson import get_ode_exprs
sympy.init_printing()
%matplotlib inline
In [ ]:
logc, logt, reduced, powsimp = False, False, False, False
odesys1 = SymbolicSys.from_callback(get_ode_exprs(logc, logt, reduced)[0],
2 if reduced else 3, 6 if reduced else 3)
odesys1.exprs
In [ ]:
odesys1._jac
Converting to sparse representation of jacobian:
In [ ]:
sparse = SymbolicSys.from_other(odesys1, sparse=True)
sparse.nnz, sparse._rowvals, sparse._colptrs
In [ ]:
sparse._jac
Converting back to dense representation of jacobian:
In [ ]:
odesys2 = SymbolicSys.from_other(sparse, sparse=False)
assert odesys2.nnz == -1
odesys2._jac
Solving the initial value problem numerically:
In [ ]:
def solve_and_plot(odesys, integrate_kw, tout=[0, 1e11], y0=[1,0,0],
params=[0.04, 1e4, 3e7], ax=None):
result = odesys.integrate(tout, y0, params, nsteps=5000,
integrator='cvode', **integrate_kw)
result.plot(ax=ax, title_info=2)
In [ ]:
import pycvodes
has_klu = pycvodes.config['KLU']
systems = [odesys1, sparse, odesys2]
kws = [{}, dict(linear_solver='klu'), {}]
fig, axes = plt.subplots(1, len(systems), figsize=(5*len(systems), 4),
subplot_kw=dict(xscale='log', yscale='log'))
for odesys, ax, kw in zip(systems, axes, kws):
if kw.get('linear_solver', '') == 'klu' and not has_klu:
continue
solve_and_plot(odesys, kw, ax=ax)