polykin.copolymerization

monomer_drift_multi

monomer_drift_multi(
    f0: FloatVectorLike,
    x: FloatVectorLike,
    r: FloatSquareMatrix,
    atol: float = 0.0001,
    rtol: float = 0.0001,
) -> FloatMatrix

Compute the monomer composition drift for a multicomponent system.

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

\[ \frac{\textup{d} f_i}{\textup{d}x} = \frac{f_i - F_i}{1 - x} \]

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

PARAMETER	DESCRIPTION
f0	Vector of initial instantaneous comonomer compositions. TYPE: FloatVectorLike(N)
x	Vector of total monomer conversion values where the drift is to be evaluated. TYPE: FloatVectorLike(M)
r	Matrix of reactivity ratios, \(r_{ij}=k_{ii}/k_{ij}\). TYPE: FloatSquareMatrix(N, N)
atol	Absolute tolerance of the ODE solver. TYPE: float DEFAULT: 0.0001
rtol	Relative tolerance of the ODE solver. TYPE: float DEFAULT: 0.0001

RETURNS	DESCRIPTION
FloatMatrix(M, N)	Matrix of monomer fraction of monomer \(i\) at the specified total monomer conversions, \(f_i(x_j)\).

See also

• inst_copolymer_multi: instantaneous copolymer composition.
• monomer_drift_multi: specific method for binary systems.

Examples:

>>> from polykin.copolymerization import monomer_drift_multi
>>> import numpy as np

Define reactivity ratio matrix.

>>> r = np.ones((3, 3))
>>> r[0, 1] = 0.2
>>> r[1, 0] = 2.3
>>> r[0, 2] = 3.0
>>> r[2, 0] = 0.9
>>> r[1, 2] = 0.4
>>> r[2, 1] = 1.5

Evaluate monomer drift.

>>> f0 = [0.5, 0.3, 0.2]
>>> x = [0.1, 0.5, 0.9, 0.99]
>>> f = monomer_drift_multi(f0, x, r)
>>> f
array([[5.19272893e-01, 2.87851432e-01, 1.92875675e-01],
       [6.38613228e-01, 2.04334321e-01, 1.57052451e-01],
       [8.31122266e-01, 5.58847454e-03, 1.63289259e-01],
       [4.98294381e-01, 1.22646553e-07, 5.01705497e-01]])

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

def monomer_drift_multi(f0: FloatVectorLike,
                        x: FloatVectorLike,
                        r: FloatSquareMatrix,
                        atol: float = 1e-4,
                        rtol: float = 1e-4
                        ) -> FloatMatrix:
    r"""Compute the monomer composition drift for a multicomponent system.

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

    $$ \frac{\textup{d} f_i}{\textup{d}x} = \frac{f_i - F_i}{1 - x} $$

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

    Parameters
    ----------
    f0 : FloatVectorLike (N)
        Vector of initial instantaneous comonomer compositions.
    x : FloatVectorLike (M)
        Vector of total monomer conversion values where the drift is to be
        evaluated.
    r : FloatSquareMatrix (N, N)
        Matrix of reactivity ratios, $r_{ij}=k_{ii}/k_{ij}$.
    atol : float
        Absolute tolerance of the ODE solver.
    rtol : float
        Relative tolerance of the ODE solver.

    Returns
    -------
    FloatMatrix (M, N)
        Matrix of monomer fraction of monomer $i$ at the specified total
        monomer conversions, $f_i(x_j)$.

    See also
    --------
    * [`inst_copolymer_multi`](inst_copolymer_multi.md): 
      instantaneous copolymer composition.
    * [`monomer_drift_multi`](monomer_drift_multi.md):
      specific method for binary systems.

    Examples
    --------
    >>> from polykin.copolymerization import monomer_drift_multi
    >>> import numpy as np

    Define reactivity ratio matrix.
    >>> r = np.ones((3, 3))
    >>> r[0, 1] = 0.2
    >>> r[1, 0] = 2.3
    >>> r[0, 2] = 3.0
    >>> r[2, 0] = 0.9
    >>> r[1, 2] = 0.4
    >>> r[2, 1] = 1.5

    Evaluate monomer drift.
    >>> f0 = [0.5, 0.3, 0.2]
    >>> x = [0.1, 0.5, 0.9, 0.99]
    >>> f = monomer_drift_multi(f0, x, r)
    >>> f
    array([[5.19272893e-01, 2.87851432e-01, 1.92875675e-01],
           [6.38613228e-01, 2.04334321e-01, 1.57052451e-01],
           [8.31122266e-01, 5.58847454e-03, 1.63289259e-01],
           [4.98294381e-01, 1.22646553e-07, 5.01705497e-01]])

    """

    N = len(f0)
    if r.shape != (N, N):
        raise ShapeError("Shape mismatch between `f0` and `r`.")

    def dfdx(x: float, f: FloatVector) -> FloatVector:
        F = inst_copolymer_multi(f=np.append(f, 1. - f.sum()), r=r)
        return (f - F[:-1]) / (1. - x + eps)

    sol = solve_ivp(dfdx,
                    t_span=(0., max(x)),
                    y0=f0[:-1],
                    t_eval=x,
                    atol=atol,
                    rtol=rtol,
                    method='LSODA')

    if sol.success:
        f = np.empty((len(x), N))
        f[:, :-1] = np.transpose(sol.y)
        f[:, -1] = 1. - np.sum(f[:, :-1], axis=1)
    else:
        raise ODESolverError(sol.message)

    return f