.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/example_asm1_ddotSO_0lvlset.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_example_asm1_ddotSO_0lvlset.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_example_asm1_ddotSO_0lvlset.py:


ASM1: equations for the 0-level set of teh 2nd time derivative of SO
============================================================================

.. GENERATED FROM PYTHON SOURCE LINES 6-43

.. code-block:: Python

    # Copyright (C) 2025 Juan Pablo Carbajal
    # Copyright (C) 2025 Mariane Yvonne Schneider
    #
    # This program is free software; you can redistribute it and/or modify
    # it under the terms of the GNU General Public License as published by
    # the Free Software Foundation; either version 3 of the License, or
    # (at your option) any later version.
    #
    # This program is distributed in the hope that it will be useful,
    # but WITHOUT ANY WARRANTY; without even the implied warranty of
    # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    # GNU General Public License for more details.
    #
    # You should have received a copy of the GNU General Public License
    # along with this program. If not, see <http://www.gnu.org/licenses/>.

    # Author: Juan Pablo Carbajal <ajuanpi@gmail.com>

    import matplotlib.pyplot as plt
    from sympy import *


    try:
        import dsafeatures
    except ModuleNotFoundError:
        import sys
        import os

        sys.path.insert(0, os.path.abspath("../.."))

    from dsafeatures.printing import *
    from dsafeatures.odemodel import ODEModel

    init_printing()

    plt.style.use('figures.mplstyle')








.. GENERATED FROM PYTHON SOURCE LINES 44-46

Select model
-------------

.. GENERATED FROM PYTHON SOURCE LINES 46-70

.. code-block:: Python

    model = ODEModel.from_csv(name='ASM1')

    X = model.states
    r_sym, r = model.rates, model.rates_expr
    M = model.coef_matrix
    O2, idxO2 = model.state_by_name("S_O2")

    # Aeration
    # Define aeration
    O2_max = Symbol(r'\text{SO}_\text{sat}', real=True, positive=True)
    max_O2_inflow = Symbol(r'\dot{q}_{\text{SO}\text{max}}', real=True, positive=True)
    aer = Function(r'\dot{q}_\text{SO}', real=True, postivive=True)(O2, O2_max, max_O2_inflow)

    # Extended system
    r_sym_ext = r_sym.col_join(Matrix([aer]))
    one_O2 = eye(len(X))[:, [idxO2]]
    M = M.row_join(one_O2)

    # First derivative
    dotX = M @ r_sym_ext
    # 2nd derivative os SO
    Jr = r_sym_ext.jacobian(X)
    ddotO2 = (M[idxO2, :] @ (Jr @ dotX))[0]








.. GENERATED FROM PYTHON SOURCE LINES 71-72

We just use numbered variables for states and paramters

.. GENERATED FROM PYTHON SOURCE LINES 72-77

.. code-block:: Python

    x = symbols(fr'x1:{len(X)+1}', real=True, positive=True)
    state_subs = dict(zip(X, x))
    p = symbols(fr'p1:{len(model.parameters())+1}', real=True, positive=True)
    param_subs = dict(zip(model.parameters(), p))








.. GENERATED FROM PYTHON SOURCE LINES 78-79

We use static symbols for the rates and their derivatives

.. GENERATED FROM PYTHON SOURCE LINES 79-84

.. code-block:: Python

    f = symbols(fr'r1:{len(r_sym_ext)+1}', real=True)
    rate_subs = dict(zip(r_sym_ext, f))
    partials = [d for d in Jr.flat() if not d.is_zero]
    df = [Symbol(fr'\,{{}}_{{{latex(state_subs[d.args[1][0]])}}}{latex(rate_subs[d.args[0]])}', real=True) for d in partials]
    dr_subs = dict(zip(partials, df))







.. GENERATED FROM PYTHON SOURCE LINES 85-86

replace these on the 2nd derivative

.. GENERATED FROM PYTHON SOURCE LINES 86-94

.. code-block:: Python

    ddO2_s = ddotO2.subs(dr_subs).subs(rate_subs).subs(param_subs)

    # collect variables until nothing changes
    ddO2_vars = [s for s in ddO2_s.free_symbols if s not in param_subs.values()]
    ddO2_s = ddO2_s.expand()
    while (expr_ := ddO2_s.collect(ddO2_vars)) != ddO2_s:
        ddO2_s = expr_








.. GENERATED FROM PYTHON SOURCE LINES 95-97

Further, we replace constant expressions with new parameter names.
To do this we need to detect all expressions that contain only parameters

.. GENERATED FROM PYTHON SOURCE LINES 97-114

.. code-block:: Python

    def has_only(expr, syms):
        """Check if expressions has only syms."""
        return not bool(expr.free_symbols - set(syms))

    def subs_expr_withall(expr, syms, x):
        if has_only(expr, syms):
            return x


    pnew = numbered_symbols(prefix='k', start=1)
    pbundle_subs = {}
    for a in preorder_traversal(ddO2_s):
        if a.is_Atom or not has_only(a, p):
            continue
        elif a not in pbundle_subs.keys():
            p_ = next(pnew)
            pbundle_subs[a] = p_







.. GENERATED FROM PYTHON SOURCE LINES 115-118

We now replace the new symbols.
However, this might leave us with repeated new constants.
Hence, we find the remaining constants and renumber them

.. GENERATED FROM PYTHON SOURCE LINES 118-124

.. code-block:: Python

    ddO2_s = ddO2_s.subs(pbundle_subs).collect(pbundle_subs.values())
    remain_ = ddO2_s.free_symbols - set(ddO2_vars)
    subs_ = {r: s for r, s in zip(remain_, numbered_symbols(prefix='C', start=1))}
    ddO2_s = ddO2_s.subs(subs_)
    pbundle_subs = {k:subs_[v] for k,v in pbundle_subs.items() if v in subs_.keys()}








.. GENERATED FROM PYTHON SOURCE LINES 125-126

Final expression

.. GENERATED FROM PYTHON SOURCE LINES 126-128

.. code-block:: Python

    ddO2_s.expand().collect(ddO2_vars).collect(pbundle_subs.values())






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    $$C_{10} \,{}_{x_{2}}r_{1} r_{2} + C_{3} \left(\,{}_{x_{6}}r_{3} r_{5} - \,{}_{x_{8}}r_{3} r_{9}\right) + C_{6} \left(\,{}_{x_{2}}r_{1} r_{7} + \,{}_{x_{5}}r_{1} r_{2} - \,{}_{x_{5}}r_{1} r_{4}\right) + \,{}_{x_{10}}r_{1} \left(C_{6} r_{6} + C_{8} \left(r_{1} + r_{2}\right) + C_{9} r_{3}\right) + \,{}_{x_{10}}r_{3} \left(C_{1} r_{3} - C_{3} r_{6} + C_{4} \left(r_{1} + r_{2}\right)\right) + \,{}_{x_{8}}r_{1} \left(C_{2} r_{1} + C_{6} r_{9} + C_{7} r_{3}\right) + \,{}_{x_{8}}r_{9} \left(- C_{3} r_{3} + C_{6} r_{1} + r_{9}\right) + r_{1} \left(C_{10} \,{}_{x_{2}}r_{1} + C_{6} \,{}_{x_{5}}r_{1} + C_{7} \,{}_{x_{8}}r_{3}\right) + r_{3} \left(- C_{3} \,{}_{x_{6}}r_{3} + C_{5} \,{}_{x_{8}}r_{3}\right)$$
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 129-129

.. code-block:: Python

    plt.show()







.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.982 seconds)

**Estimated memory usage:**  88 MB


.. _sphx_glr_download_auto_examples_example_asm1_ddotSO_0lvlset.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: example_asm1_ddotSO_0lvlset.ipynb <example_asm1_ddotSO_0lvlset.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: example_asm1_ddotSO_0lvlset.py <example_asm1_ddotSO_0lvlset.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: example_asm1_ddotSO_0lvlset.zip <example_asm1_ddotSO_0lvlset.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_