<jupyter_start><jupyter_text>One-dimensional Heat Equation<jupyter_code>import numpy as np
import scipy.linalg as la
import scipy.sparse as sparse
import matplotlib.pyplot as plt

import rom_operator_inference as opinf
# Matplotlib customizations.
plt.rc("figure", dpi=300, figsize=(9,3))
plt.rc("font", family="serif")
np.random.seed(10)<jupyter_output><empty_output><jupyter_text>This example is based on the Operator Inference problem for the 1D heat equation described in {cite}`PW2016OperatorInference`.  Problem Statement Let $\Omega = [0,L]\subset \mathbb{R}$ be the spatial domain indicated by the variable $x$, and let $[0,T]\subset\mathbb{R}$ be the time domain with variable $t$.We consider the one-dimensional heat equation with non-homogeneous Dirichlet boundary conditions,\begin{align*}    \frac{\partial}{\partial t} q(x,t) - \frac{\partial^2}{\partial x^2}q(x,t) &= 0    & x &\in\Omega,\quad t\in[0,T],    \\    q(0,t) = q(L,t) &= u(t)    & t &\in[0,T],    \\    q(x,0) = \big(e^{\alpha(x - 1)} + e^{-\alpha x} &- e^{-\alpha}\big)u(0)    & x &\in \Omega.\end{align*} Let $\{x_i\}_{i=0}^{n+1}$ be an equidistant grid of $n+2$ points on $\Omega$, i.e.,\begin{align*}    0 &= x_0 < x_1 < \cdots < x_n < x_{n+1} = L    &    &\text{and}    &    \delta x &= \frac{L}{n+1} = x_{i+1} - x_{i},\quad i=1,\ldots,n-1.\end{align*}Since the boundary conditions prescribe $q(x_0,t) = q(x_{n+1},t) = 1$, we wish to compute the state vector $\mathbf{q}(t) = \begin{bmatrix} q(x_1,t) & \cdots & q(x_n,t)\end{bmatrix}^{\top}\in\mathbb{R}^n$ for various $t\in[0,T]$. Introducing the finite difference approximation\begin{align*}    \frac{\partial^2}{\partial x^2}q(x,t) &\approx \frac{q(x-\delta x,t) - 2q(x,t) + q(x+\delta x,t)}{(\delta x)^2}    % &    % \Longrightarrow&    % &    % \frac{\partial^2}{\partial x^2}q(x_i,t) &\approx \frac{q(x_{i-1},t) - 2q(x_{i},t) + q(x_{i+1},t)}{(\delta x)^2}    % \\    &    &\Longrightarrow    &    \frac{\partial^2}{\partial x^2}q_{i} &\approx \frac{q_{i-1} - 2q_{i} + q_{i+1}}{(\delta x)^2},\end{align*}we obtain the semi-discrete linear system$$\frac{\text{d}}{\text{d}t}\mathbf{q}(t) = \mathbf{A}\mathbf{q}(t) + \mathbf{B}u(t),$$where\begin{align*}    \mathbf{A} &= \frac{1}{(\delta x)^2}\left[\begin{array}{ccccc}        -2 & 1 & & & \\        1 & -2 & 1 & & \\        & \ddots & \ddots & \ddots & \\        & & 1 & -2 & 1 \\        & & & 1 & -2 \\    \end{array}\right] \in\mathbb{R}^{n\times n},    &    \mathbf{B} &= \frac{1}{(\delta x)^2}\left[\begin{array}{c}        1 \\ 0 \\ \vdots \\ 0 \\ 1    \end{array}\right]\in\mathbb{R}^{n}.\end{align*}  Snapshot Data Generation For simplicity, let $L = T = 1$ and $u(t) = 1$.We begin by simulating the full-order system described above with a maximal time step size $\delta t = 10^{-3}$, resulting in $k = 10^3+1$ time steps (1000 steps past the initial condition).The result is the snapshot matrix $\mathbf{Q}\in\mathbb{R}^{n\times k}$, where the $j$th column is the solution trajectory at time $t_j$.We also compute the time derivative at each snapshot, obtaining $\dot{\mathbf{Q}}\in\mathbb{R}^{n\times k}$.<jupyter_code># Construct the spatial domain.
L = 1                           # Spatial domain length.
n = 2**7 - 1                    # Spatial grid size.
x_all = np.linspace(0, L, n+2)  # Full spatial grid.
x = x_all[1:-1]                 # Interior spatial grid (where q is unknown).
dx = x[1] - x[0]                # Spatial resolution.

# Construct the temporal domain.
T = 1                           # Temporal domain length (final simulation time).
k = T*10**3 + 1                 # Temporal grid size.
t = np.linspace(0, T, k)        # Temporal grid.
dt = t[1] - t[0]                # Temporal resolution.

print(f"Spatial step size δx = {dx}")
print(f"Temporal step size δt = {dt}")
# Construct state matrix A.
dx2inv = 1 / dx**2
diags = np.array([1, -2, 1]) * dx2inv
A = sparse.diags(diags, [-1,0,1], (n,n))

# Construct input matrix B.
B = np.zeros_like(x)
B[0], B[-1] = dx2inv, dx2inv
u = lambda t: np.ones_like(t)   # Input function u(t) = 1.
U = u(t)                        # Inputs over the time domain.

# Construct the initial condition.
alpha = 100
q0 = np.exp(alpha*(x-1)) + np.exp(-alpha*x) - np.exp(-alpha)

print(f"shape of A:\t{A.shape}")
print(f"shape of B:\t{B.shape}")
print(f"shape of q0:\t{q0.shape}")<jupyter_output><empty_output><jupyter_text>Since this is a diffusive problem, we will use the Implicit (Backward) Euler method for solving the ODEs.For the problem $\frac{\text{d}}{\text{d}t}\mathbf{q}(t) = \mathbf{f}(t, \mathbf{q}(t), \mathbf{u}(t))$, the method is defined by$$    \mathbf{q}_{j+1} = \mathbf{q}_{j} + \delta t \mathbf{f}(t_{j+1},\mathbf{q}_{j+1},u_{j+1}).$$With the form $\mathbf{f}(t,\mathbf{q}(t),u(t)) = \mathbf{A}\mathbf{q}(t) + \mathbf{B}u(t)$, this becomes$$    \mathbf{q}_{j+1} = (I - \delta t \mathbf{A})^{-1}\left(\mathbf{q}_{j} + \delta t \mathbf{B} u_{j+1}\right).$$<jupyter_code>def implicit_euler(t, q0, A, B, U):
    """Solve the system

        dq / dt = Aq(t) + Bu(t),    q(0) = q0,

    over a uniform time domain via the Implicit Euler method.

    Parameters
    ----------
    t : (k,) ndarray
        Uniform time array over which to solve the ODE.
    q0 : (n,) ndarray
        Initial condition.
    A : (n,n) ndarray
        State matrix.
    B : (n,) or (n,1) ndarray
        Input matrix.
    U : (k,) ndarray
        Inputs over the time array.

    Returns
    -------
    q : (n,k) ndarray
        Solution to the ODE at time t; that is, q[:,j] is the
        computed solution corresponding to time t[j].
    """
    # Check and store dimensions.
    k = len(t)
    n = len(q0)
    B = np.ravel(B)
    assert A.shape == (n,n)
    assert B.shape == (n,)
    assert U.shape == (k,)
    I = np.eye(n)

    # Check that the time step is uniform.
    dt = t[1] - t[0]
    assert np.allclose(np.diff(t), dt)

    # Factor I - dt*A for quick solving at each time step.
    factored = la.lu_factor(I - dt*A)

    # Solve the problem at each time step.
    q = np.empty((n,k))
    q[:,0] = q0.copy()
    for j in range(1,k):
        q[:,j] = la.lu_solve(factored, q[:,j-1] + dt*B*U[j])

    return q
# Compute snapshots by solving the equation with implicit_euler().
Q = implicit_euler(t, q0, A, B, U)

# Also compute time derivatives (dq/dt) at each snapshot.
Qdot = A @ Q + B.reshape((-1,1))*U

print(f"shape of Q:\t{Q.shape}")
print(f"shape of Qdot:\t{Qdot.shape}")<jupyter_output><empty_output><jupyter_text>Finally, we visualize the snapshots to get a sense of how the solution looks qualitatively.<jupyter_code>def plot_heat_data(Z, title="Snapshot Data"):
    fig, [ax1,ax2] = plt.subplots(1, 2, figsize=(12,4))

    # Plot a few snapshots.
    color = iter(plt.cm.viridis(np.linspace(.25, 1, 5)))
    for j in [0, 20, 80, 160, 640]:
        q_all = np.concatenate([[1], Z[:,j], [1]])  # Pad with boundary conditions.
        ax1.plot(x_all, q_all, color=next(color), label=f"q(x,t_{j})")
    ax1.set_xlim(0, 1)
    ax1.set_xlabel("x")
    ax1.set_ylabel("q(x,t)")
    ax1.legend(loc="lower right", fontsize=8, bbox_to_anchor=(1.01,.05))

    # Plot all snapshots in space and time.
    xx, tt = np.meshgrid(x, t, indexing="ij")
    cdata = ax2.pcolormesh(xx, tt, Z, shading="nearest", cmap="magma")
    plt.colorbar(cdata, ax=ax2, extend="both")
    ax2.set_xlabel("space")
    ax2.set_ylabel("time")

    fig.suptitle(title)
    plt.show()
plot_heat_data(Q)<jupyter_output><empty_output><jupyter_text>Reduced Model Construction Now that we have snapshot data $\mathbf{Q}$, we can construct a POD basis $\mathbf{V}_r$ to use in the construction of the ROM.There are a few ways to make an informed choice of $r$; in this example, we examine the _relative projection error_, defined by$$\text{err}_\text{projection} = \frac{||\mathbf{Q} - \mathbf{V}_r \mathbf{V}_r^{\top}\mathbf{Q}||_F}{||\mathbf{Q}||_F}.$$<jupyter_code># Compute the largest possible basis and all singular values.
V, svdvals = opinf.pre.pod_basis(Q)

# For 1 ≤ r < 20, calculate the projection error of the snapshot matrix.
rs = np.arange(1, 21)
projection_errors = [opinf.pre.projection_error(Q, V[:,:r])[1] for r in rs]

# Plot the errors.
plt.semilogy(rs, projection_errors, "C0d-")
plt.axhline(1e-5, color='k', lw=.5, alpha=.5)
plt.axvline(8, color='k', lw=.5, alpha=.5)
plt.xlim(rs[0], rs[-1])
plt.xticks(rs[::2])
plt.xlabel("Reduced dimension")
plt.ylabel("Relative projection error")
plt.show()<jupyter_output><empty_output><jupyter_text>We choose $r=8$, the smallest number of basis vectors for which the projection error is less than $10^{-5}$.<jupyter_code># Print absolute and relative projection errors for r = 8.
Vr = V[:,:8]
opinf.pre.projection_error(Q, Vr)<jupyter_output><empty_output><jupyter_text>Now we can learn the reduced model with Operator Inference.Because the full-order model is of the form $\frac{\text{d}}{\text{d}t}\mathbf{q}(t) = \mathbf{A}\mathbf{q}(t) + \mathbf{B}u(t)$, we specify a reduced model form of $\frac{\text{d}}{\text{d}t}\widehat{\mathbf{q}}(t) = \hat{\mathbf{A}}\widehat{\mathbf{q}}(t) + \hat{\mathbf{B}}u(t)$ (`modelform="AB"`).<jupyter_code># Train and run the model.
inferred_rom = opinf.ContinuousOpInfROM(modelform="AB")
inferred_rom.fit(Vr, Q, Qdot, U)
q0_ = Vr.T @ q0                                     # Project the initial condition.
Q_ROM_inferred = Vr @ implicit_euler(t, q0_, inferred_rom.A_.entries, inferred_rom.B_.entries, U)<jupyter_output><empty_output><jupyter_text>For the error analysis, since we used the projection error to determine $r$, we use a similar measure to evaluate the state error:$$\text{err}_\text{state} = \frac{||\mathbf{Q} - \mathbf{Q}_\text{ROM}||_F}{||\mathbf{Q}||_F}.$$<jupyter_code>opinf.post.frobenius_error(Q, Q_ROM_inferred)[1]
plot_heat_data(Q_ROM_inferred, "Output of Inferred ROM")<jupyter_output><empty_output><jupyter_text>We can also check how well we did relative to the projection error.<jupyter_code>relative_projection_error = opinf.post.lp_error(Q, Vr @ Vr.T @ Q, normalize=True)[1]

def plot_relative_errors_over_time(Zlist, labels):
    colors = ["C0", "C3"]
    plt.semilogy(t, relative_projection_error, "C1", label="Projection Error")
    for Z,label,c in zip(Zlist, labels, colors[:len(Zlist)]):
        relative_l2_error = opinf.post.lp_error(Q, Z, normalize=True)[1]
        plt.semilogy(t, relative_l2_error, c, label=label)

    plt.xlabel("t")
    plt.ylabel("Relative Error")
    plt.legend(loc="upper right", edgecolor="none")
    plt.show()
plot_relative_errors_over_time([Q_ROM_inferred], ["OpInf ROM Error"])<jupyter_output><empty_output><jupyter_text>Comparison with Intrusive Methods Under some idealized assumptions, the operators learned through Operator Inference "converge" in a sense to the corresponding operators obtained through intrusive projection,\begin{align*}    \widetilde{\mathbf{A}} &= \mathbf{V}_{r}^{\top} \mathbf{A} \mathbf{V}_{r}^{\top},    &    \widetilde{\mathbf{B}} &= \mathbf{V}_{r}^{\top}\mathbf{B}.\end{align*}We construct the (intrusive) reduced model corresponding to these projected reduced operators for comparison.<jupyter_code>intrusive_model = opinf.ContinuousOpInfROM("AB")
intrusive_model.fit(Vr, None, None, known_operators={"A":A, "B":B})
Q_ROM_intrusive = Vr @ implicit_euler(t, q0_, intrusive_model.A_.entries, intrusive_model.B_.entries, U)
plot_heat_data(Q_ROM_intrusive, "Output of Intrusive ROM")
plot_relative_errors_over_time([Q_ROM_inferred, Q_ROM_intrusive],
                               ["OpInf ROM Error", "Intrusive ROM Error"])
la.norm(Q_ROM_intrusive - Q_ROM_inferred) / la.norm(Q_ROM_intrusive)<jupyter_output><empty_output><jupyter_text>In this case, **the inferred and projected ROMs give essentially the same result**.However, the inferred ROM successfully emulates the FOM **without explicit knowledge of the operators** $\mathbf{A}$ **and** $\mathbf{B}$.Before moving forward, let's see how the dimension $r$ affects the accuracy of the ROM.<jupyter_code>V, vals = opinf.pre.pod_basis(Q)

def run_trial(r):
    Vr = V[:,:r]

    # Construct and simulate the intrusive ROM.
    intrusive_model = opinf.ContinuousOpInfROM("AB").fit(Vr, None, None,
                                                         known_operators={"A":A, "B":B})
    Q_ROM_intrusive = Vr @ implicit_euler(t, Vr.T @ q0,
                                          intrusive_model.A_.entries, intrusive_model.B_.entries, U)

    # Construct and simulate the inferred ROM.
    inferred_rom = opinf.ContinuousOpInfROM("AB").fit(Vr, Q, Qdot, U)
    Q_ROM_inferred = Vr @ implicit_euler(t, Vr.T @ q0,
                                         inferred_rom.A_.entries, inferred_rom.B_.entries, U)

    # Calculate errors.
    projection_error = opinf.pre.projection_error(Q, Vr)[1]
    intrusive_error = opinf.post.frobenius_error(Q, Q_ROM_intrusive)[1]
    inference_error = opinf.post.frobenius_error(Q, Q_ROM_inferred)[1]

    return projection_error, intrusive_error, inference_error
def plot_state_error(rmax):
    rs = np.arange(1, rmax+1)
    err_projection, err_intrusive, err_inference = zip(*[run_trial(r) for r in rs])

    plt.semilogy(rs, err_projection, 'C1-', label="projection error")
    plt.semilogy(rs, err_intrusive, 'C3+-', label="intrusive", mew=2)
    plt.semilogy(rs, err_inference, 'C0o-', label="inference-based", mfc='none', mec='C0', mew=1.5)

    plt.xlim(rs.min(), rs.max())
    plt.xticks(rs, [str(int(r)) for r in rs])
    plt.xlabel("Reduced dimension")
    plt.ylabel("Relative error")
    plt.legend(loc="upper right", fontsize=14, framealpha=1)
    plt.grid(ls=':')
    plt.show()
plot_state_error(14)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17032806092560046, 'fraction_numerical': 0.017867603983596953, 'mean_word_length': 3.6484002722940776, 'pattern_counts': {'":': 4, '<': 50, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '450689', 'license': 'MIT', 'n_tokens_mistral': 5142, 'n_tokens_neox': 4938, 'n_words': 1238, 'path': 'docs/content/examples/heat_1D.ipynb', 'repo_name': 'mtezzele/rom-operator-inference-Python3', 'size': '13656'}