polykin.copolymerization

monomer_drift_binary

monomer_drift_binary(
    f10: Union[float, FloatVectorLike],
    x: FloatVectorLike,
    r1: float,
    r2: float,
    atol: float = 0.0001,
    rtol: float = 0.0001,
    method: Literal["LSODA", "RK45"] = "LSODA",
) -> FloatMatrix

Compute the monomer composition drift for a binary system.

In a closed binary system, the drift in monomer composition is given by the solution of the following differential equation:

\[ \frac{\textup{d} f_1}{\textup{d}x} = \frac{f_1 - F_1}{1 - x} \]

with initial condition \(f_1(0)=f_{1,0}\), where \(f_1\) and \(F_1\) are, respectively, the instantaneous comonomer and copolymer composition of M1, and \(x\) is the total molar monomer conversion.

PARAMETER	DESCRIPTION
f10	Initial molar fraction of M1, \(f_{1,0}=f_1(0)\). TYPE: float | FloatVectorLike(N)
x	Total monomer conversion values where the drift is to be evaluated. TYPE: FloatVectorLike(M)
r1	Reactivity ratio of M1. TYPE: float | FloatArray
r2	Reactivity ratio of M2. TYPE: float | FloatArray
atol	Absolute tolerance of ODE solver. TYPE: float DEFAULT: 0.0001
rtol	Relative tolerance of ODE solver. TYPE: float DEFAULT: 0.0001
method	ODE solver. TYPE: Literal['LSODA', 'RK45'] DEFAULT: 'LSODA'

RETURNS	DESCRIPTION
FloatMatrix(M, N)	Monomer fraction of M1 at a given conversion, \(f_1(x)\).

See also

• monomer_drift_multi: generic method for multicomponent systems.

Examples:

>>> from polykin.copolymerization import monomer_drift_binary

An example with f10 as scalar.

>>> f1 = monomer_drift_binary(f10=0.5, x=[0.1, 0.5, 0.9], r1=0.16, r2=0.70)
>>> f1
array([0.51026241, 0.57810678, 0.87768138])

An example with f10 as list.

>>> f1 = monomer_drift_binary(f10=[0.2, 0.8], x=[0.1, 0.5, 0.9],
...                           r1=0.16, r2=0.70)
>>> f1
array([[0.19841009, 0.18898084, 0.15854395],
       [0.82315475, 0.94379024, 0.99996457]])

Source code in src/polykin/copolymerization/binary.py

def monomer_drift_binary(f10: Union[float, FloatVectorLike],
                         x: FloatVectorLike,
                         r1: float,
                         r2: float,
                         atol: float = 1e-4,
                         rtol: float = 1e-4,
                         method: Literal['LSODA', 'RK45'] = 'LSODA'
                         ) -> FloatMatrix:
    r"""Compute the monomer composition drift for a binary system.

    In a closed binary system, the drift in monomer composition is given by
    the solution of the following differential equation:

    $$ \frac{\textup{d} f_1}{\textup{d}x} = \frac{f_1 - F_1}{1 - x} $$

    with initial condition $f_1(0)=f_{1,0}$, where $f_1$ and $F_1$ are,
    respectively, the instantaneous comonomer and copolymer composition of
    M1, and $x$ is the total molar monomer conversion.

    Parameters
    ----------
    f10 : float | FloatVectorLike (N)
        Initial molar fraction of M1, $f_{1,0}=f_1(0)$.
    x : FloatVectorLike (M)
        Total monomer conversion values where the drift is to be evaluated.
    r1 : float | FloatArray
        Reactivity ratio of M1.
    r2 : float | FloatArray
        Reactivity ratio of M2.
    atol : float
        Absolute tolerance of ODE solver.
    rtol : float
        Relative tolerance of ODE solver.
    method : Literal['LSODA', 'RK45']
        ODE solver.

    Returns
    -------
    FloatMatrix (M, N)
        Monomer fraction of M1 at a given conversion, $f_1(x)$.

    See also
    --------
    * [`monomer_drift_multi`](monomer_drift_multi.md):
      generic method for multicomponent systems.

    Examples
    --------
    >>> from polykin.copolymerization import monomer_drift_binary

    An example with f10 as scalar.
    >>> f1 = monomer_drift_binary(f10=0.5, x=[0.1, 0.5, 0.9], r1=0.16, r2=0.70)
    >>> f1
    array([0.51026241, 0.57810678, 0.87768138])

    An example with f10 as list.
    >>> f1 = monomer_drift_binary(f10=[0.2, 0.8], x=[0.1, 0.5, 0.9],
    ...                           r1=0.16, r2=0.70)
    >>> f1
    array([[0.19841009, 0.18898084, 0.15854395],
           [0.82315475, 0.94379024, 0.99996457]])
    """

    if isinstance(f10, (int, float)):
        f10 = [f10]

    sol = solve_ivp(df1dx,
                    (0., max(x)),
                    f10,
                    args=(r1, r2),
                    t_eval=x,
                    method=method,
                    vectorized=True,
                    atol=atol,
                    rtol=rtol)

    if sol.success:
        result = sol.y
        result = np.maximum(0., result)
        if result.shape[0] == 1:
            result = result[0]
    else:
        raise ODESolverError(sol.message)

    return result