In Chapters 8–11 you learned that many lumped systems can be written as \( \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t) \). For distributed-parameter systems, the “state” is typically a field \( x(\xi,t) \) over space \( \xi\in(0,L) \). A mathematically faithful way to retain the state-space viewpoint is to embed the field into a Hilbert space such as \( H=L^2(0,L) \) with inner product \( \langle f,g\rangle = \int_0^L f(\xi)g(\xi)\,d\xi \).

We then write an abstract evolution equation \( x(t)\in H \):

\[ \dot{x}(t)=A x(t) + B u(t),\qquad y(t)=C x(t),\qquad x(0)=x_0. \]

Here \( A \) is typically a differential operator (e.g., Laplacian for diffusion), and \( B \), \( C \) encode how actuators and sensors couple into the PDE. The central subtlety: \( A \) is often unbounded, so solutions must be defined with functional-analytic tools (semigroups) rather than naive matrix exponentials.

flowchart TD
  S["Physical PDE model (field state x(xi,t))"] --> H["Choose Hilbert space H (e.g., L2(0,L))"]
  H --> Aop["Define operator A (spatial derivatives + BCs)"]
  Aop --> IO["Define input/output operators B, C (actuators/sensors)"]
  IO --> SG["Semigroup solution x(t)=T(t)x0 + int T(t-s)Bu(s) ds"]
  SG --> DIS["Discretize (FD/FEM) -> large ODE xdot = A_N x + B_N u"]
  DIS --> CTRL["Design controller/observer (e.g., LQR/Kalman) on reduced model"]
  CTRL --> VAL["Validate on refined discretization / full PDE"]
        

Let \( A:D(A)\subset H \to H \) be a densely-defined linear operator. If \( A \) generates a strongly continuous (C₀) semigroup \( (T(t))_{t\ge 0} \), then the homogeneous system \( \dot{x}(t)=Ax(t) \) has solution \( x(t)=T(t)x_0 \). With inputs, the appropriate notion is the mild solution:

\[ x(t)=T(t)x_0 + \int_0^t T(t-s)\,B u(s)\,ds. \]

Key theorem (existence/uniqueness, mild form). If \( A \) generates a C₀ semigroup and \( u\in L^2(0,T;U) \) (input space \( U \)) and \( B\in\mathcal{L}(U,H) \) (bounded), then for every \( x_0\in H \) there exists a unique mild solution \( x\in C([0,T];H) \).

Proof sketch. Use semigroup boundedness \( \|T(t)\|\le M e^{\omega t} \) and estimate:

\[ \left\| \int_0^t T(t-s)Bu(s)\,ds \right\| \le \int_0^t \|T(t-s)\|\,\|B\|\,\|u(s)\|\,ds \le M\|B\|\int_0^t e^{\omega(t-s)}\|u(s)\|\,ds, \]

which shows the integral term is well-defined in \( H \). Uniqueness follows from linearity and semigroup properties. When \( B \) is unbounded (typical for boundary control), one replaces boundedness by admissibility (Section 3).

Generator conditions (intuition). For many PDEs, \( A \) is dissipative and maximal, so by the Lumer–Phillips theorem it generates a contraction semigroup \( \|T(t)\|\le 1 \).

Two canonical couplings:

• Distributed actuation: an input density \( b(\xi)u(t) \) enters the PDE inside the domain, often giving a bounded operator \( B:U\to H \).
• Boundary actuation: the control enters through boundary conditions (Dirichlet/Neumann/Robin). In the abstract state-space, this typically yields an unbounded control operator. Well-posedness then uses admissible control operators.

A standard 1D diffusion example with distributed actuation:

\[ \frac{\partial x}{\partial t}(\xi,t)=\alpha\frac{\partial^2 x}{\partial \xi^2}(\xi,t)+b(\xi)u(t), \qquad x(0,t)=0,\;x(L,t)=0. \]

If \( b\in L^2(0,L) \), then \( B: \mathbb{R}\to L^2(0,L) \) defined by \( (Bu)(\xi)=b(\xi)u \) is bounded. For boundary actuation (e.g., \( x(0,t)=u(t) \)), \( B \) becomes unbounded in \( L^2 \), and admissibility replaces boundedness:

\[ \exists T>0,\;\exists c_T>0:\quad \int_0^T \left\| T(T-s)B u(s)\right\|_H^2\,ds \le c_T \int_0^T \|u(s)\|_U^2\,ds. \]

Similarly, pointwise sensing (e.g., \( y(t)=x(\xi_s,t) \)) is not bounded on \( L^2(0,L) \) and is handled via admissible observation operators.

flowchart TD
  U["Input u(t)"] --> BOP["Actuator operator B"]
  BOP --> DYN["Infinite-dim dynamics: xdot = A x + B u"]
  DYN --> COP["Sensor operator C"]
  COP --> Y["Output y(t)"]
  Y --> FB["Controller: u(t) = -K x_hat(t)"]
  FB --> U
        

PDE stability is often proved via energy functionals, the infinite-dimensional analog of quadratic Lyapunov functions \( V(\mathbf{x})=\mathbf{x}^\top P\mathbf{x} \).

For the diffusion equation with Dirichlet boundaries, define \( E(t)=\tfrac{1}{2}\|x(t)\|_{L^2}^2 = \tfrac{1}{2}\int_0^L x(\xi,t)^2\,d\xi \). Using integration by parts (and boundary terms vanishing):

\[ \frac{dE}{dt} = \int_0^L x\,x_t\,d\xi = \alpha\int_0^L x\,x_{\xi\xi}\,d\xi = -\alpha\int_0^L (x_\xi)^2\,d\xi \le 0. \]

Thus the semigroup is contractive in \( L^2 \) and solutions decay (in an appropriate sense). For wave-like PDEs the undamped energy is conserved; stabilization requires damping (internal or boundary).

Operator viewpoint. The inequality above corresponds to \( \langle Ax,x\rangle \le 0 \) (dissipativity). Maximal dissipativity implies \( A \) generates a contraction semigroup.

In Lesson 2 you studied separation of variables and eigenfunction expansions. Write \( x(\xi,t)=\sum_{n=1}^\infty x_n(t)\phi_n(\xi) \), where \( \{\phi_n\} \) are eigenfunctions of the spatial operator \( \mathcal{A} \). For 1D diffusion with Dirichlet boundaries:

\[ \phi_n(\xi)=\sqrt{\frac{2}{L} }\sin\!\left(\frac{n\pi\xi}{L}\right), \qquad \lambda_n=-\alpha\left(\frac{n\pi}{L}\right)^2. \]

Projecting the PDE yields countably many ODEs:

\[ \dot{x}_n(t)=\lambda_n x_n(t) + b_n u(t), \qquad b_n=\langle b,\phi_n\rangle. \]

Truncation. Approximating by the first \( N \) modes gives \( \dot{\mathbf{x} }_N = \mathbf{A}_N\mathbf{x}_N+\mathbf{B}_N u \). For diffusion, high-frequency modes decay fast because \( \lambda_n \to -\infty \), which often makes low-order truncations accurate for control. A simple bound (informal but instructive) for the tail energy is:

\[ \left\| (I-P_N) x(t) \right\|_{L^2}^2 = \sum_{n>N} x_n(t)^2 \le e^{2\lambda_{N+1} t}\sum_{n>N} x_n(0)^2 \le e^{2\lambda_{N+1} t}\|x(0)\|_{L^2}^2, \]

where \( P_N \) is the projection onto the first \( N \) modes. This highlights a crucial spatio-temporal control idea: diffusion “filters” spatial high frequencies.

The finite-dimensional concepts extend using operator-valued Gramians. For bounded \( B \) and a semigroup \( T(t) \), define the reachability (controllability) Gramian over \( [0,T] \):

\[ W_T = \int_0^T T(t)\,B B^\ast\,T(t)^\ast\,dt, \]

where \( B^\ast \) is the adjoint and \( T(t)^\ast \) is the adjoint semigroup. Exact controllability in infinite dimensions is delicate and often fails for diffusion (compactness issues), while approximate controllability may hold under mild conditions on \( b_n \).

Observability similarly uses \( W_T^{obs}=\int_0^T T(t)^\ast C^\ast C T(t)\,dt \). Duality between controllability and observability remains central, but proofs require functional analysis (adjoint operators, dense ranges, compact embeddings).

The simplest link is state feedback. In the abstract form, one seeks \( u(t)=-Kx(t) \) so that \( A_{cl}=A-BK \) generates an exponentially stable semigroup. For bounded \( K \) this parallels finite-dimensional design.

A canonical optimal-control formulation (infinite-dimensional LQR) is:

\[ \min_{u(\cdot)}\; J(u)=\int_0^\infty \left(\langle Qx(t),x(t)\rangle_H + \langle Ru(t),u(t)\rangle_U\right)dt, \quad \dot{x}=Ax+Bu. \]

Under standard stabilizability/detectability assumptions, the value functional is quadratic \( V(x)=\langle Px,x\rangle \) with a self-adjoint positive operator \( P \) solving an operator Riccati equation:

\[ A^\ast P + P A - P B R^{-1} B^\ast P + Q = 0, \qquad u^\ast(t) = -R^{-1}B^\ast P x(t). \]

This equation is the infinite-dimensional analog of the matrix CARE you will use after discretization. In practice, most engineering workflows compute an approximate controller by (i) discretizing the PDE into a large ODE, (ii) solving a matrix Riccati equation, (iii) validating the closed-loop on refined meshes.

If only \( y(t)=Cx(t) \) is measured, a Luenberger-type observer in Hilbert spaces is:

\[ \dot{\hat{x} }(t)=A\hat{x}(t)+Bu(t)+L\big(y(t)-C\hat{x}(t)\big), \]

where \( L \) is chosen so that \( A-LC \) is stable (in semigroup sense). For diffusion-like systems, sensor placement (which spatial points/regions are measured) strongly affects observability because some modes may be weakly seen (small \( \langle c,\phi_n\rangle \)).

A standard bridge is the method of lines: discretize space (FD/FEM), keep time continuous, obtain \( \dot{\mathbf{x} }=\mathbf{A}_N\mathbf{x}+\mathbf{B}_N\mathbf{u} \). This enables reuse of finite-dimensional tools (stability, LQR, observers) with caution:

• Discretization must respect the PDE’s energy structure (e.g., dissipativity) to avoid artificial instability.
• Boundary control/sensing may require special discretizations because the true operators are unbounded.
• Controller validity should be checked under mesh refinement (convergence of closed-loop behavior).

We implement a 1D heat equation with distributed actuation, discretized via finite differences. Then we design a finite-dimensional LQR controller on the discretized model. This is a practical approximation of the infinite-dimensional LQR idea in Section 7.

Chapter19_Lesson5.py

# Chapter19_Lesson5.py
# Spatio-Temporal Dynamics: 1D heat equation discretization + (finite-dimensional) LQR

import numpy as np
from scipy.linalg import solve_continuous_are

def laplacian_1d_dirichlet(N: int, L: float):
    """Second-derivative matrix on N interior points with Dirichlet boundaries."""
    dx = L / (N + 1)
    main = -2.0 * np.ones(N)
    off  =  1.0 * np.ones(N - 1)
    D2 = (np.diag(main) + np.diag(off, 1) + np.diag(off, -1)) / (dx * dx)
    x = dx * (np.arange(1, N + 1))  # interior grid
    return D2, x, dx

def gaussian_shape(x, x0=0.25, sigma=0.08):
    return np.exp(-0.5*((x-x0)/sigma)**2)

def simulate(A, B, K, x0, dt=5e-3, T=5.0):
    steps = int(T/dt)
    x = x0.copy()
    hist = np.zeros((steps+1, x.size))
    E = np.zeros(steps+1)
    hist[0] = x
    E[0] = 0.5 * float(x.T @ x)
    for k in range(steps):
        u = -float(K @ x)  # scalar input
        xdot = A @ x + B.flatten() * u
        x = x + dt * xdot
        hist[k+1] = x
        E[k+1] = 0.5 * float(x.T @ x)
    return hist, E

def main():
    # PDE: x_t = alpha x_xx + b(x) u(t),  x(0,t)=x(L,t)=0
    alpha = 0.12
    L = 1.0
    N = 60

    D2, grid, dx = laplacian_1d_dirichlet(N, L)
    A = alpha * D2

    # Distributed actuator shape b(x) (bounded input operator in the PDE idealization)
    b = gaussian_shape(grid, x0=0.25, sigma=0.07)
    b = b / np.linalg.norm(b)  # normalize
    B = b.reshape(-1, 1)

    # Finite-dimensional LQR: minimize \int (x^T Q x + u^T R u) dt
    Q = np.eye(N)
    R = np.array([[2e-3]])

    # Continuous-time algebraic Riccati equation: A^T P + P A - P B R^{-1} B^T P + Q = 0
    P = solve_continuous_are(A, B, Q, R)
    K = np.linalg.solve(R, B.T @ P)  # 1xN

    # Initial condition: a smooth bump
    x0 = np.exp(-80.0*(grid-0.75)**2)

    hist, E = simulate(A, B, K, x0, dt=2e-3, T=4.0)

    print("Grid points N =", N, "dx =", dx)
    print("Closed-loop gain K shape:", K.shape)
    print("Energy E(0)   =", E[0])
    print("Energy E(end) =", E[-1])
    print("Energy decay factor =", E[-1]/E[0])

if __name__ == "__main__":
    main()

Chapter19_Lesson5.cpp

// Chapter19_Lesson5.cpp
// Discretized 1D heat equation + discrete-time Riccati iteration (DARE) using Eigen
//
// Build with (example):
//   g++ -O2 -std=c++17 Chapter19_Lesson5.cpp -I /path/to/eigen -o heat_lqr

#include <iostream>
#include <cmath>
#include <Eigen/Dense>

using Eigen::MatrixXd;
using Eigen::VectorXd;

static MatrixXd laplacian_dirichlet(int N, double L) {
    double dx = L / (N + 1.0);
    MatrixXd D2 = MatrixXd::Zero(N, N);
    for (int i = 0; i < N; ++i) {
        D2(i, i) = -2.0;
        if (i - 1 >= 0) D2(i, i - 1) = 1.0;
        if (i + 1 <  N) D2(i, i + 1) = 1.0;
    }
    return D2 / (dx * dx);
}

static VectorXd grid_interior(int N, double L) {
    double dx = L / (N + 1.0);
    VectorXd x(N);
    for (int i = 0; i < N; ++i) x(i) = dx * (i + 1.0);
    return x;
}

static VectorXd gaussian_shape(const VectorXd& x, double x0, double sigma) {
    VectorXd b(x.size());
    for (int i = 0; i < x.size(); ++i) {
        double z = (x(i) - x0) / sigma;
        b(i) = std::exp(-0.5 * z * z);
    }
    return b;
}

// Simple DARE fixed-point iteration:
//   P_{k+1} = Q + A^T P A - A^T P B (R + B^T P B)^{-1} B^T P A
static MatrixXd dare_iterate(const MatrixXd& A, const MatrixXd& B,
                             const MatrixXd& Q, const MatrixXd& R,
                             int iters = 5000, double tol = 1e-10) {
    MatrixXd P = Q;
    for (int k = 0; k < iters; ++k) {
        MatrixXd BtPB = B.transpose() * P * B;
        MatrixXd S = R + BtPB;
        MatrixXd K = S.inverse() * (B.transpose() * P * A); // (m x n)
        MatrixXd Pnext = Q + A.transpose() * P * A - A.transpose() * P * B * K;
        double err = (Pnext - P).norm() / (P.norm() + 1e-12);
        P = Pnext;
        if (err < tol) break;
    }
    return P;
}

int main() {
    // Semi-discrete model: x_dot = A x + B u  -> Euler: x_{k+1} = Ad x_k + Bd u_k
    const double alpha = 0.12;
    const double L = 1.0;
    const int N = 60;

    MatrixXd D2 = laplacian_dirichlet(N, L);
    MatrixXd A = alpha * D2;

    VectorXd xgrid = grid_interior(N, L);
    VectorXd b = gaussian_shape(xgrid, 0.25, 0.07);
    b /= b.norm();
    MatrixXd B = b;                 // (N x 1)

    // Discretize
    const double dt = 2e-3;
    MatrixXd Ad = MatrixXd::Identity(N, N) + dt * A;
    MatrixXd Bd = dt * B;

    MatrixXd Q = MatrixXd::Identity(N, N);
    MatrixXd R(1,1); R(0,0) = 2e-3;

    MatrixXd P = dare_iterate(Ad, Bd, Q, R, 20000, 1e-12);

    MatrixXd S = R + Bd.transpose() * P * Bd;
    MatrixXd K = S.inverse() * (Bd.transpose() * P * Ad); // (1 x N)

    // Initial condition
    VectorXd x0(N);
    for (int i = 0; i < N; ++i) {
        double z = xgrid(i) - 0.75;
        x0(i) = std::exp(-80.0 * z * z);
    }

    // Simulate
    int steps = int(4.0 / dt);
    VectorXd x = x0;
    double E0 = 0.5 * x.squaredNorm();
    for (int k = 0; k < steps; ++k) {
        double u = -(K * x)(0,0);
        x = Ad * x + Bd * u;
    }
    double E1 = 0.5 * x.squaredNorm();

    std::cout << "N=" << N << " dt=" << dt << "\n";
    std::cout << "Energy E(0)=" << E0 << "  E(end)=" << E1
              << "  decay=" << (E1/(E0+1e-18)) << "\n";
    return 0;
}

Chapter19_Lesson5.java

// Chapter19_Lesson5.java
// Discretized 1D heat equation + discrete-time Riccati iteration (DARE) using EJML
//
// Maven dependency (example):
//   org.ejml:ejml-simple:0.43

import org.ejml.simple.SimpleMatrix;

public class Chapter19_Lesson5 {

    static SimpleMatrix laplacianDirichlet(int N, double L) {
        double dx = L / (N + 1.0);
        SimpleMatrix D2 = new SimpleMatrix(N, N);
        for (int i = 0; i < N; i++) {
            D2.set(i, i, -2.0);
            if (i - 1 >= 0) D2.set(i, i - 1, 1.0);
            if (i + 1 < N)  D2.set(i, i + 1, 1.0);
        }
        return D2.divide(dx * dx);
    }

    static SimpleMatrix gridInterior(int N, double L) {
        double dx = L / (N + 1.0);
        SimpleMatrix x = new SimpleMatrix(N, 1);
        for (int i = 0; i < N; i++) x.set(i, 0, dx * (i + 1.0));
        return x;
    }

    static SimpleMatrix gaussianShape(SimpleMatrix x, double x0, double sigma) {
        int N = x.numRows();
        SimpleMatrix b = new SimpleMatrix(N, 1);
        for (int i = 0; i < N; i++) {
            double z = (x.get(i, 0) - x0) / sigma;
            b.set(i, 0, Math.exp(-0.5 * z * z));
        }
        return b;
    }

    static SimpleMatrix dareIterate(SimpleMatrix A, SimpleMatrix B,
                                   SimpleMatrix Q, SimpleMatrix R,
                                   int iters, double tol) {
        SimpleMatrix P = Q.copy();
        for (int k = 0; k < iters; k++) {
            SimpleMatrix BtPB = B.transpose().mult(P).mult(B);
            SimpleMatrix S = R.plus(BtPB);
            SimpleMatrix K = S.invert().mult(B.transpose().mult(P).mult(A)); // (m x n)
            SimpleMatrix Pnext = Q.plus(A.transpose().mult(P).mult(A))
                                  .minus(A.transpose().mult(P).mult(B).mult(K));
            double err = Pnext.minus(P).normF() / (P.normF() + 1e-12);
            P = Pnext;
            if (err < tol) break;
        }
        return P;
    }

    public static void main(String[] args) {
        final double alpha = 0.12;
        final double L = 1.0;
        final int N = 60;

        SimpleMatrix D2 = laplacianDirichlet(N, L);
        SimpleMatrix A = D2.scale(alpha);

        SimpleMatrix xgrid = gridInterior(N, L);
        SimpleMatrix b = gaussianShape(xgrid, 0.25, 0.07);
        b = b.divide(b.normF());
        SimpleMatrix B = b; // (N x 1)

        final double dt = 2e-3;
        SimpleMatrix Ad = SimpleMatrix.identity(N).plus(A.scale(dt));
        SimpleMatrix Bd = B.scale(dt);

        SimpleMatrix Q = SimpleMatrix.identity(N);
        SimpleMatrix R = new SimpleMatrix(1, 1);
        R.set(0, 0, 2e-3);

        SimpleMatrix P = dareIterate(Ad, Bd, Q, R, 20000, 1e-12);
        SimpleMatrix S = R.plus(Bd.transpose().mult(P).mult(Bd));
        SimpleMatrix K = S.invert().mult(Bd.transpose().mult(P).mult(Ad)); // (1 x N)

        // Initial condition
        SimpleMatrix x = new SimpleMatrix(N, 1);
        for (int i = 0; i < N; i++) {
            double z = xgrid.get(i, 0) - 0.75;
            x.set(i, 0, Math.exp(-80.0 * z * z));
        }

        int steps = (int)(4.0 / dt);
        double E0 = 0.5 * x.dot(x);
        for (int k = 0; k < steps; k++) {
            double u = -K.mult(x).get(0, 0);
            x = Ad.mult(x).plus(Bd.scale(u));
        }
        double E1 = 0.5 * x.dot(x);

        System.out.println("N=" + N + " dt=" + dt);
        System.out.println("Energy E(0)=" + E0 + "  E(end)=" + E1 + "  decay=" + (E1/(E0+1e-18)));
    }
}

Chapter19_Lesson5.m

% Chapter19_Lesson5.m
% Spatio-Temporal Dynamics: method of lines for 1D heat equation + LQR feedback
%
% PDE: x_t = alpha x_xx + b(x) u(t),  x(0,t)=x(L,t)=0.
% Discretize to x_dot = A x + B u, then use LQR: u = -K x.

clear; clc;

alpha = 0.12;
L = 1.0;
N = 60;                 % interior nodes
dx = L/(N+1);

% Laplacian (Dirichlet)
e = ones(N,1);
D2 = spdiags([e -2*e e], [-1 0 1], N, N) / (dx^2);
A = alpha * full(D2);

xgrid = dx*(1:N)';

% Actuator shape (Gaussian)
x0 = 0.25; sigma = 0.07;
b = exp(-0.5*((xgrid-x0)/sigma).^2);
b = b / norm(b);
B = b; % (N x 1)

Q = eye(N);
R = 2e-3;

% Continuous-time LQR
% Requires Control System Toolbox: [K,~,~] = lqr(A,B,Q,R);
[K,~,~] = lqr(A,B,Q,R);

% Initial condition: smooth bump
x_init = exp(-80*(xgrid-0.75).^2);

% Simulate with forward Euler (for teaching purposes)
dt = 2e-3;
T  = 4.0;
steps = floor(T/dt);

x = x_init;
E = zeros(steps+1,1);
E(1) = 0.5*(x'*x);

for k = 1:steps
    u = -K*x;
    xdot = A*x + B*u;
    x = x + dt*xdot;
    E(k+1) = 0.5*(x'*x);
end

fprintf('Energy E(0)   = %.6e\n', E(1));
fprintf('Energy E(end) = %.6e\n', E(end));
fprintf('Decay factor  = %.6e\n', E(end)/E(1));

% Optional plot
figure; plot(linspace(0,T,steps+1), E);
xlabel('t'); ylabel('E(t)=0.5*||x||^2'); grid on;
title('Closed-loop energy decay (finite-dimensional LQR on PDE discretization)');

Chapter19_Lesson5.nb

(* Chapter19_Lesson5.nb *)
(* Spatio-Temporal Dynamics: 1D heat equation discretization + discrete Riccati feedback *)

ClearAll[alpha, L, N, dx, D2, A, xgrid, b, B, dt, Ad, Bd, Q, R, P, K];

alpha = 0.12;
L = 1.0;
N = 60;
dx = L/(N + 1);

(* Laplacian with Dirichlet boundaries on interior nodes *)
D2 = SparseArray[
   {Band[{1, 1}] -> (-2), Band[{1, 2}] -> 1, Band[{2, 1}] -> 1},
   {N, N}
] / dx^2;

A = alpha Normal[D2];
xgrid = dx Range[N];

(* Actuator shape *)
b = Exp[-1/2 ((xgrid - 0.25)/0.07)^2];
b = b/Norm[b];
B = Transpose[{b}]; (* N x 1 *)

(* Discretize *)
dt = 2.*10^-3;
Ad = IdentityMatrix[N] + dt A;
Bd = dt B;

Q = IdentityMatrix[N];
R = { {2.*10^-3} };

(* Discrete algebraic Riccati equation solver *)
P = DiscreteRiccatiSolve[{Ad, Bd}, Q, R];
K = Inverse[R + Transpose[Bd].P.Bd].Transpose[Bd].P.Ad; (* 1 x N *)

(* Initial condition *)
x = Exp[-80 (xgrid - 0.75)^2];

E0 = 0.5 x.x;

steps = Floor[4.0/dt];
Do[
  u = - (K.x)[[1]];
  x = Ad.x + Flatten[Bd]*u;
, {k, 1, steps}];

E1 = 0.5 x.x;

Print["Energy E(0)=", E0];
Print["Energy E(end)=", E1];
Print["Decay factor=", E1/(E0 + 10^-18)];

Problem 1 (Mild solution and semigroup property): Assume \( A \) generates a C₀ semigroup \( T(t) \) on a Hilbert space \( H \) and \( B\in\mathcal{L}(U,H) \). Show that the function \( x(t)=T(t)x_0+\int_0^t T(t-s)Bu(s)\,ds \) satisfies the variation-of-constants identity and is unique.

Solution: Define \( \Phi(t)=\int_0^t T(t-s)Bu(s)\,ds \). Using the semigroup property \( T(t)T(\tau)=T(t+\tau) \) and a change of variables,

\[ x(t+h)=T(h)x(t)+\int_t^{t+h} T(t+h-s)Bu(s)\,ds. \]

This is the integral form of the dynamics. Uniqueness follows by taking two solutions, subtracting, and using Grönwall-type arguments in the semigroup norm bound (or linearity plus density of \( D(A) \)).

Problem 2 (Energy decay for diffusion): For \( x_t=\alpha x_{\xi\xi} \) on \( (0,L) \) with \( x(0,t)=x(L,t)=0 \), prove \( \frac{d}{dt}\frac{1}{2}\int_0^L x^2\,d\xi \le 0 \).

Solution: Multiply by \( x \) and integrate:

\[ \frac{d}{dt}\frac{1}{2}\int_0^L x^2\,d\xi = \int_0^L x x_t\,d\xi = \alpha\int_0^L x x_{\xi\xi}\,d\xi = \alpha\left[x x_\xi\right]_0^L - \alpha\int_0^L (x_\xi)^2\,d\xi = -\alpha\int_0^L (x_\xi)^2\,d\xi \le 0, \]

since \( x(0,t)=x(L,t)=0 \Rightarrow [x x_\xi]_0^L=0 \).

Problem 3 (Mode dynamics with distributed actuation): Let \( x(\xi,t)=\sum_{n\ge1} x_n(t)\phi_n(\xi) \) be an eigenfunction expansion of the diffusion operator with Dirichlet boundaries. For input \( b(\xi)u(t) \), derive the scalar mode ODEs and an expression for \( b_n \).

Solution: Project onto \( \phi_n \): \( \dot{x}_n=\langle x_t,\phi_n\rangle \). Using \( \mathcal{A}\phi_n=\lambda_n\phi_n \) and orthonormality:

\[ \dot{x}_n(t)=\lambda_n x_n(t)+\langle b,\phi_n\rangle u(t), \qquad b_n=\langle b,\phi_n\rangle = \int_0^L b(\xi)\phi_n(\xi)\,d\xi. \]

Problem 4 (Discrete Riccati iteration): Consider \( x_{k+1}=A_d x_k + B_d u_k \) and cost \( \sum_{k=0}^\infty (x_k^\top Q x_k + u_k^\top R u_k) \), with \( Q\succeq 0 \), \( R\succ 0 \). Show that if \( P \) solves the DARE, the optimal gain is \( K=(R+B_d^\top P B_d)^{-1}B_d^\top P A_d \).

Solution: Using dynamic programming with value function \( V(x)=x^\top P x \), minimize \( x^\top Q x + u^\top R u + (A_d x + B_d u)^\top P (A_d x + B_d u) \) over \( u \). Setting derivative w.r.t. \( u \) to zero yields \( (R+B_d^\top P B_d)u + B_d^\top P A_d x = 0 \), hence the stated \( K \).

Problem 5 (Truncation tail bound for diffusion): For diffusion with eigenvalues \( \lambda_n<0 \), show that the tail energy satisfies \( \sum_{n>N} x_n(t)^2 \le e^{2\lambda_{N+1}t}\sum_{n>N} x_n(0)^2 \).

Solution: Each unforced mode obeys \( \dot{x}_n=\lambda_n x_n \Rightarrow x_n(t)=e^{\lambda_n t}x_n(0) \). Since \( \lambda_n \le \lambda_{N+1} \) for all \( n>N \),

\[ \sum_{n>N} x_n(t)^2 = \sum_{n>N} e^{2\lambda_n t} x_n(0)^2 \le e^{2\lambda_{N+1}t}\sum_{n>N} x_n(0)^2. \]

We formulated spatio-temporal PDE dynamics as infinite-dimensional state-space systems on Hilbert spaces, introduced semigroup well-posedness and mild solutions, and discussed how actuation/measurement can be bounded (distributed) or unbounded (boundary) and handled via admissibility. We connected modal decompositions to control-relevant truncation and showed how matrix Riccati-based feedback on PDE discretizations provides a practical bridge to distributed control theory.