In Chapter 19 Lesson 1–2, we introduced PDE models and modal decompositions (separation of variables). Numerical discretization instead replaces the spatial field by finitely many degrees of freedom. After discretization, the PDE becomes a (typically large) ODE/DAE that fits the state-space framework of Chapters 8–11.

Consider a scalar PDE on \( x\in[0,L] \): \( u_t = \mathcal{A}u + s(x,t) \), with boundary conditions. Spatial discretization produces a vector state \( \mathbf{u}(t)\in\mathbb{R}^n \) and an operator approximation \( \mathbf{A}_h \) such that

\[ \dot{\mathbf{u} }(t) = \mathbf{A}_h\,\mathbf{u}(t) + \mathbf{b}_h(t), \quad \mathbf{u}(0)=\mathbf{u}_0. \]

\( h \): spatial mesh size, \( n \): number of spatial degrees of freedom (often \( n\sim 1/h \)). The time-discretization advances this high-order ODE using explicit or implicit integrators, producing a fully discrete recursion \( \mathbf{u}^{k+1}=\Phi(\mathbf{u}^k) \).

flowchart TD
  A["PDE model on x in [0,L]"] --> B["Choose spatial grid/mesh (size h)"]
  B --> C["Approximate spatial operator: u_x, u_xx, ..."]
  C --> D["Obtain MOL ODE: u_dot = A_h*u + b_h(t)"]
  D --> E["Choose time integrator: explicit/implicit"]
  E --> F["Simulate & analyze stability/accuracy"]
  F --> G["Use as high-order state-space model for control/estimation"]
        

Two key points (often confused in practice):

• Physical stability (PDE well-posedness/dissipation) is about the true continuous model.
• Numerical stability is about the discrete update map \( \Phi \) and depends on \( h \), \( \Delta t \), and the scheme.

Let a uniform grid be \( x_i = i h \): for \( i=0,1,\dots,N \), where \( h=L/N \). Denote \( u_i(t)=u(x_i,t) \). FD methods replace derivatives by algebraic combinations (stencils) of nearby samples.

2.1 First derivative

Forward, backward, and central differences:

\[ u_x(x_i,t) \approx \frac{u_{i+1}(t)-u_i(t)}{h} \quad (\text{forward}),\qquad u_x(x_i,t) \approx \frac{u_i(t)-u_{i-1}(t)}{h} \quad (\text{backward}), \\ u_x(x_i,t) \approx \frac{u_{i+1}(t)-u_{i-1}(t)}{2h} \quad (\text{central}). \]

2.2 Second derivative

A standard second-order approximation of \( u_{xx} \) is the three-point stencil:

\[ u_{xx}(x_i,t) \approx \frac{u_{i+1}(t)-2u_i(t)+u_{i-1}(t)}{h^2}. \]

2.3 Truncation error via Taylor expansion (proof sketch)

Expand about \( x_i \):

\[ u(x_i \pm h,t)=u(x_i,t) \pm h u_x(x_i,t) + \frac{h^2}{2}u_{xx}(x_i,t) \pm \frac{h^3}{6}u_{xxx}(x_i,t) + \frac{h^4}{24}u_{xxxx}(x_i,t) + \mathcal{O}(h^5). \]

Subtracting the “+” and “−” expansions eliminates even powers, giving for the central first derivative:

\[ \frac{u(x_i+h,t)-u(x_i-h,t)}{2h} = u_x(x_i,t) + \frac{h^2}{6}u_{xxx}(x_i,t) + \mathcal{O}(h^4), \]

hence second-order accuracy. Adding the expansions eliminates odd powers and yields:

\[ \frac{u(x_i+h,t)-2u(x_i,t)+u(x_i-h,t)}{h^2} = u_{xx}(x_i,t) + \frac{h^2}{12}u_{xxxx}(x_i,t) + \mathcal{O}(h^4), \]

so the second-derivative stencil is also second-order.

2.4 Discrete operators as matrices

On interior nodes with homogeneous Dirichlet boundary conditions, the FD Laplacian yields a tridiagonal matrix (for interior size \( n=N-1 \)):

\[ \mathbf{A}_h = \frac{1}{h^2} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & 1 & -2 & 1 \\ 0 & \cdots & 0 & 1 & -2 \end{bmatrix}. \]

For Dirichlet BCs, \( \mathbf{A}_h \) is symmetric negative definite. Its eigenvalues are explicitly known:

\[ \lambda_j(\mathbf{A}_h)= -\frac{4}{h^2}\sin^2\!\left(\frac{j\pi}{2(n+1)}\right),\quad j=1,\dots,n. \]

In particular, the most negative eigenvalue satisfies \( |\lambda_{\max}| \approx 4/h^2 \), which directly controls explicit time-step restrictions for diffusion.

3.1 Heat equation as a large ODE

For the 1D heat equation \( u_t = \alpha u_{xx} \) with Dirichlet BCs, FD in space yields the MOL system:

\[ \dot{\mathbf{u} }(t) = \alpha \mathbf{A}_h\mathbf{u}(t) + \mathbf{b}_h(t), \]

where \( \mathbf{b}_h(t) \) encodes boundary and source contributions. This is precisely a high-dimensional linear state-space model.

3.2 Theta-method family (explicit & implicit)

The theta method for \( \dot{\mathbf{u} }=\mathbf{f}(t,\mathbf{u}) \) is:

\[ \mathbf{u}^{k+1} = \mathbf{u}^k + \Delta t\left[\theta\,\mathbf{f}(t_{k+1},\mathbf{u}^{k+1}) + (1-\theta)\,\mathbf{f}(t_k,\mathbf{u}^k)\right], \quad \theta\in[0,1]. \]

For linear diffusion MOL, this becomes a linear solve each step:

\[ \left(\mathbf{I}-\theta\Delta t\,\alpha\mathbf{A}_h\right)\mathbf{u}^{k+1} = \left(\mathbf{I}+(1-\theta)\Delta t\,\alpha\mathbf{A}_h\right)\mathbf{u}^{k} + \Delta t\,\mathbf{b}_{\theta}. \]

Special cases: \( \theta=0 \) (forward Euler), \( \theta=1/2 \) (Crank–Nicolson), \( \theta=1 \) (backward Euler).

3.3 Stability from spectra (matrix viewpoint)

Suppose boundary/source terms are absent so \( \dot{\mathbf{u} }=\alpha\mathbf{A}_h\mathbf{u} \). For forward Euler:

\[ \mathbf{u}^{k+1} = \left(\mathbf{I}+\Delta t\,\alpha\mathbf{A}_h\right)\mathbf{u}^k. \]

Stability requires the spectral radius \( \rho(\mathbf{I}+\Delta t\,\alpha\mathbf{A}_h) \le 1 \). Since \( \mathbf{A}_h \) is diagonalizable with real negative eigenvalues, this reduces to requiring \( |1+\Delta t\,\alpha\lambda_j|\le 1 \) for all \( j \), i.e.

\[ -2 \le \Delta t\,\alpha\lambda_j \le 0 \quad \Rightarrow\quad \Delta t \le \frac{2}{\alpha|\lambda_{\max}|}. \]

Using \( |\lambda_{\max}|\approx 4/h^2 \) yields the classical restriction:

\[ \Delta t \le \frac{h^2}{2\alpha}. \]

3.4 Stability: von Neumann viewpoint and CFL-type restrictions

For explicit FD schemes on uniform grids, Fourier-mode (von Neumann) analysis gives the same result as the spectral argument above. For FTCS heat, define \( r=\alpha\Delta t/h^2 \); a sufficient stability condition is:

\[ r = \frac{\alpha\Delta t}{h^2} \le \frac{1}{2}. \]

Implicit choices (e.g., \( \theta\ge 1/2 \)) are A-stable for linear diffusion operators: stable for any \( \Delta t \) on the test equation \( y'=\lambda y \) with \( \Re(\lambda) < 0 \). (Accuracy still depends on \( \Delta t \).)

3.5 Consistency and the stability → convergence principle

A cornerstone theorem for linear IVPs is Lax–Richtmyer: for a properly posed linear initial-value problem, a consistent FD scheme converges if and only if it is stable. Conceptually:

\[ \text{Consistency} + \text{Stability} \Rightarrow \text{Convergence}. \]

Practically: (i) compute truncation error (consistency), (ii) check stability (spectral, energy, or Fourier), then (iii) verify convergence by mesh refinement.

4.1 Weak form for Poisson (integration by parts)

For \( -u_{xx}=f \) on \( (0,L) \) with Dirichlet BC, multiply by test \( v \) (vanishing at the boundary) and integrate by parts:

\[ \int_0^L u'(x)\,v'(x)\,dx = \int_0^L f(x)\,v(x)\,dx, \quad \forall v \in H_0^1(0,L). \]

4.2 Galerkin approximation

Choose a mesh \( 0=x_0 < x_1 < \dots < x_N=L \) and a finite-dimensional subspace \( V_h \subset H_0^1(0,L) \) spanned by hat functions \( \{\varphi_i\} \). Approximate \( u_h(x)=\sum_i U_i \varphi_i(x) \). Testing with \( v=\varphi_j \) yields:

\[ \sum_i U_i \underbrace{\int_0^L \varphi_i'(x)\varphi_j'(x)\,dx}_{K_{ji} } = \underbrace{\int_0^L f(x)\varphi_j(x)\,dx}_{F_j} \quad \Rightarrow\quad \mathbf{K}\mathbf{U}=\mathbf{F}. \]

4.3 Element matrices (local-to-global assembly)

On an element of length \( h_e \), the P1 stiffness matrix is:

\[ \mathbf{K}_e = \frac{1}{h_e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}. \]

Assembly sums these local contributions into the global sparse matrix. Dirichlet BCs are enforced by elimination or penalty, yielding an SPD system.

4.4 FEM energy norm error bound (Céa-type estimate)

Let \( a(u,v)=\int_0^L u'v'\,dx \) and \( \|w\|_a^2=a(w,w) \). Under standard assumptions (coercivity and continuity), the Galerkin solution \( u_h\in V_h \) satisfies:

\[ \|u-u_h\|_a \le \frac{C_{\text{cont} } }{C_{\text{coer} } }\,\inf_{v_h\in V_h}\|u-v_h\|_a. \]

This says FEM is quasi-optimal: it performs (up to a constant) as well as the best approximation in the chosen space. For P1 elements on a uniform mesh, one obtains (for sufficiently smooth solutions) typical rates: \( \|u-u_h\|_{H^1}=\mathcal{O}(h) \), \( \|u-u_h\|_{L^2}=\mathcal{O}(h^2) \).

4.5 Transient diffusion: mass matrix and semi-discrete ODE

For the heat equation weak form:

\[ \int_0^L u_t\,v\,dx + \alpha\int_0^L u'\,v'\,dx = \int_0^L s\,v\,dx, \]

FEM in space yields:

\[ \mathbf{M}\dot{\mathbf{U} }(t) + \alpha\mathbf{K}\mathbf{U}(t) = \mathbf{F}(t), \quad \mathbf{M}_{ji}=\int_0^L \varphi_i\varphi_j\,dx. \]

flowchart TD
  A["Start from PDE"] --> B["Multiply by test function v and integrate"]
  B --> C["Integrate by parts -> weak form"]
  C --> D["Choose mesh + basis (hat functions)"]
  D --> E["Assemble element matrices -> global K and M"]
  E --> F["Apply boundary conditions"]
  F --> G["Semi-discrete system: M*U_dot + K*U = F"]
  G --> H["Time stepping on matrix ODE/DAE"]
        

• Geometry: FD is simplest on structured grids; FEM naturally supports irregular meshes and complex domains.
• Derivatives vs energies: FD approximates derivatives pointwise; FEM approximates weak/energy forms, often improving boundary/heterogeneity handling.
• Linear algebra: both yield sparse systems; diffusion yields SPD matrices after BC enforcement (FD Laplacian; FEM stiffness/mass).
• Time-step restrictions: explicit diffusion typically needs \( \Delta t=\mathcal{O}(h^2) \); implicit removes stiffness at the cost of solves.
• Control perspective: discretization yields high-dimensional state-space models; discretization error becomes “unmodeled dynamics” for controllers/observers.

File: Chapter19_Lesson3.py

"""
Chapter19_Lesson3.py
Finite Difference (FD) and Finite Element (FE) Approximations (Conceptual Level)

This script provides:
  1) 1D heat equation discretization via finite differences (method of lines),
     with explicit FTCS and implicit theta-method (e.g., Crank–Nicolson).
  2) 1D Poisson / steady diffusion via linear finite elements (Galerkin FEM).

Dependencies:
  - numpy
  - scipy (sparse matrices and solvers)
  - matplotlib (optional, for plots)

Author: (course material generator)
"""

from __future__ import annotations

import numpy as np
from dataclasses import dataclass
from typing import Callable, Tuple

from scipy.sparse import diags, csc_matrix, identity
from scipy.sparse.linalg import spsolve


@dataclass
class Grid1D:
    L: float = 1.0
    N: int = 101  # number of nodes including boundaries

    def __post_init__(self) -> None:
        if self.N < 3:
            raise ValueError("N must be >= 3.")
        self.x = np.linspace(0.0, self.L, self.N)
        self.dx = self.x[1] - self.x[0]


def laplacian_dirichlet_1d(N: int, dx: float) -> csc_matrix:
    """
    FD Laplacian on interior nodes with homogeneous Dirichlet boundaries.
    Unknowns are u[1:-1] (size N-2).
    """
    n = N - 2
    main = -2.0 * np.ones(n) / (dx * dx)
    off = 1.0 * np.ones(n - 1) / (dx * dx)
    A = diags([off, main, off], offsets=[-1, 0, 1], format="csc")
    return A


def heat_fd_explicit_ftcs(
    u0: np.ndarray,
    alpha: float,
    grid: Grid1D,
    dt: float,
    steps: int,
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Solve u_t = alpha u_xx with FD in space and forward Euler in time (FTCS).
    Dirichlet boundaries are taken from u0[0] and u0[-1] and kept fixed.

    Stability (1D): r = alpha*dt/dx^2 <= 1/2  (classical sufficient condition).
    """
    u = u0.copy()
    dx = grid.dx
    r = alpha * dt / (dx * dx)
    if r > 0.5:
        print(f"[warning] FTCS may be unstable: r = {r:.4f} > 0.5")

    us = [u.copy()]
    ts = [0.0]
    for n in range(steps):
        un = u.copy()
        # interior update
        u[1:-1] = un[1:-1] + r * (un[2:] - 2.0 * un[1:-1] + un[:-2])
        us.append(u.copy())
        ts.append((n + 1) * dt)
    return np.array(ts), np.array(us)


def heat_fd_theta_method(
    u0: np.ndarray,
    alpha: float,
    grid: Grid1D,
    dt: float,
    steps: int,
    theta: float = 0.5,
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Theta-method on the method-of-lines system:
        u_t = alpha * A u + b(t), with Dirichlet BC.

    In FD form for interior vector v = u[1:-1]:
        (I - theta*dt*alpha*A) v^{n+1} = (I + (1-theta)*dt*alpha*A) v^n + dt*rhs_bc

    theta = 0   -> forward Euler  (explicit)
    theta = 1/2 -> Crank–Nicolson (A-stable, 2nd order in time)
    theta = 1   -> backward Euler (A-stable, 1st order in time)
    """
    u = u0.copy()
    dx = grid.dx
    N = grid.N
    A = laplacian_dirichlet_1d(N, dx)

    n = N - 2
    I = identity(n, format="csc")
    LHS = (I - theta * dt * alpha * A)
    RHS = (I + (1.0 - theta) * dt * alpha * A)

    # fixed Dirichlet values
    uL = float(u0[0])
    uR = float(u0[-1])

    us = [u.copy()]
    ts = [0.0]

    for k in range(steps):
        v = u[1:-1].copy()

        # Boundary contribution for Laplacian: A*v assumes uL=uR=0,
        # so we add the missing terms.
        bc = np.zeros(n)
        bc[0] += uL / (dx * dx)
        bc[-1] += uR / (dx * dx)

        rhs = RHS @ v + dt * alpha * bc
        v_new = spsolve(LHS, rhs)
        u[1:-1] = v_new
        us.append(u.copy())
        ts.append((k + 1) * dt)

    return np.array(ts), np.array(us)


def fem_poisson_1d_linear(
    Nel: int,
    L: float,
    f: Callable[[np.ndarray], np.ndarray],
    u0: float = 0.0,
    uL: float = 0.0,
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Solve -u'' = f(x) on (0,L) with Dirichlet BC u(0)=u0, u(L)=uL
    using linear (P1) FEM on a uniform mesh with Nel elements.

    Weak form:
        find u in H^1_0 + lift such that ∫ u' v' dx = ∫ f v dx  for all v in H^1_0.

    Assembly uses element stiffness:
        Ke = (1/h) [[ 1, -1],
                    [-1,  1]]
    and load with midpoint quadrature on each element.
    """
    if Nel < 1:
        raise ValueError("Nel must be >= 1.")
    x = np.linspace(0.0, L, Nel + 1)
    h = x[1] - x[0]
    Nn = Nel + 1  # nodes

    K = np.zeros((Nn, Nn))
    F = np.zeros(Nn)

    for e in range(Nel):
        i, j = e, e + 1
        Ke = (1.0 / h) * np.array([[1.0, -1.0], [-1.0, 1.0]])
        xm = 0.5 * (x[i] + x[j])
        fe = f(np.array([xm]))[0]
        # consistent load for linear basis with midpoint quadrature:
        # ∫ f phi_i dx ≈ f(xm) * h/2 for each local node
        Fe = fe * (h / 2.0) * np.array([1.0, 1.0])

        K[i:i+2, i:i+2] += Ke
        F[i:i+2] += Fe

    # Dirichlet BC via elimination (simple for teaching)
    fixed = {0: u0, Nn - 1: uL}
    free = np.array([k for k in range(Nn) if k not in fixed])

    # modify RHS for known values
    for idx, val in fixed.items():
        F = F - K[:, idx] * val

    Kff = K[np.ix_(free, free)]
    Ff = F[free]
    uf = np.linalg.solve(Kff, Ff)

    u = np.zeros(Nn)
    u[free] = uf
    for idx, val in fixed.items():
        u[idx] = val
    return x, u


def main() -> None:
    # --- Heat equation example ---
    grid = Grid1D(L=1.0, N=101)
    alpha = 0.05

    # initial condition with Dirichlet boundaries fixed to 0
    u0 = np.sin(np.pi * grid.x)
    u0[0] = 0.0
    u0[-1] = 0.0

    # explicit (stable) and Crank–Nicolson
    dt_exp = 0.4 * (grid.dx * grid.dx) / alpha  # r=0.4 <= 0.5
    steps = 200

    t_exp, u_exp = heat_fd_explicit_ftcs(u0, alpha, grid, dt_exp, steps)
    t_cn, u_cn = heat_fd_theta_method(u0, alpha, grid, dt_exp, steps, theta=0.5)

    print("Heat FD done. Final time:", t_exp[-1])

    # --- Poisson FEM example ---
    f = lambda x: (np.pi ** 2) * np.sin(np.pi * x)  # so exact u = sin(pi x)
    x_fem, u_fem = fem_poisson_1d_linear(Nel=40, L=1.0, f=f, u0=0.0, uL=0.0)
    u_exact = np.sin(np.pi * x_fem)

    err_inf = np.max(np.abs(u_fem - u_exact))
    print("FEM Poisson max error:", err_inf)

    # Optional plot (comment out if running headless)
    try:
        import matplotlib.pyplot as plt

        plt.figure()
        plt.plot(grid.x, u_exp[-1], label="FD explicit (final)")
        plt.plot(grid.x, u_cn[-1], "--", label="FD Crank–Nicolson (final)")
        plt.legend()
        plt.xlabel("x")
        plt.ylabel("u")
        plt.title("1D Heat Equation: FD in Space")
        plt.show()

        plt.figure()
        plt.plot(x_fem, u_fem, "o-", label="FEM P1")
        plt.plot(x_fem, u_exact, "--", label="exact")
        plt.legend()
        plt.xlabel("x")
        plt.ylabel("u")
        plt.title("1D Poisson: Linear FEM")
        plt.show()
    except Exception as e:
        print("Plot skipped:", e)


if __name__ == "__main__":
    main()
      

File: Chapter19_Lesson3.cpp

/*
Chapter19_Lesson3.cpp
Finite Difference (FD) and Finite Element (FE) Approximations (Conceptual Level)

This C++ example implements:
  1) 1D heat equation u_t = alpha u_xx via FD in space and theta-method in time.
  2) 1D Poisson -u'' = f(x) via linear FEM (Galerkin) on uniform mesh.

Dependencies (recommended):
  - Eigen (Dense + Sparse)
    https://eigen.tuxfamily.org/

Compile (example):
  g++ -O2 -std=c++17 Chapter19_Lesson3.cpp -I /path/to/eigen -o lesson19_3

Author: (course material generator)
*/

#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>

#include <Eigen/Dense>
#include <Eigen/Sparse>

using Eigen::VectorXd;
using Eigen::MatrixXd;
using Eigen::SparseMatrix;
using Eigen::Triplet;

struct Grid1D {
    double L;
    int N; // nodes including boundaries
    std::vector<double> x;
    double dx;

    Grid1D(double L_, int N_) : L(L_), N(N_) {
        x.resize(N);
        for (int i = 0; i < N; ++i) x[i] = (L * i) / (N - 1);
        dx = x[1] - x[0];
    }
};

SparseMatrix<double> laplacian_dirichlet_1d(int N, double dx) {
    // interior unknowns: size n = N-2
    int n = N - 2;
    std::vector<Triplet<double>> T;
    T.reserve(3 * n);

    double s = 1.0 / (dx * dx);
    for (int i = 0; i < n; ++i) {
        T.emplace_back(i, i, -2.0 * s);
        if (i - 1 >= 0) T.emplace_back(i, i - 1, 1.0 * s);
        if (i + 1 < n)  T.emplace_back(i, i + 1, 1.0 * s);
    }
    SparseMatrix<double> A(n, n);
    A.setFromTriplets(T.begin(), T.end());
    return A;
}

void heat_fd_theta(
    const VectorXd& u0_full,
    double alpha,
    const Grid1D& grid,
    double dt,
    int steps,
    double theta,
    std::vector<double>& t_out,
    std::vector<VectorXd>& u_out
) {
    int N = grid.N;
    int n = N - 2;

    SparseMatrix<double> A = laplacian_dirichlet_1d(N, grid.dx);
    SparseMatrix<double> I(n, n);
    I.setIdentity();

    SparseMatrix<double> LHS = (I - theta * dt * alpha * A);
    SparseMatrix<double> RHS = (I + (1.0 - theta) * dt * alpha * A);

    // Dirichlet boundaries fixed
    double uL = u0_full(0);
    double uR = u0_full(N - 1);

    // interior vector
    VectorXd v = u0_full.segment(1, n);

    // Pre-factorize LHS
    Eigen::SparseLU<SparseMatrix<double>> solver;
    solver.analyzePattern(LHS);
    solver.factorize(LHS);

    t_out.clear(); u_out.clear();
    t_out.reserve(steps + 1);
    u_out.reserve(steps + 1);

    auto pack_full = [&](const VectorXd& vin) {
        VectorXd u(N);
        u.setZero();
        u(0) = uL;
        u(N - 1) = uR;
        u.segment(1, n) = vin;
        return u;
    };

    t_out.push_back(0.0);
    u_out.push_back(pack_full(v));

    VectorXd bc(n); bc.setZero();
    double s = 1.0 / (grid.dx * grid.dx);
    bc(0)     += uL * s;
    bc(n - 1) += uR * s;

    for (int k = 0; k < steps; ++k) {
        VectorXd rhs = RHS * v + dt * alpha * bc;
        v = solver.solve(rhs);

        t_out.push_back((k + 1) * dt);
        u_out.push_back(pack_full(v));
    }
}

std::pair<std::vector<double>, VectorXd> fem_poisson_p1(
    int Nel, double L,
    double (*f)(double),
    double u0, double uL
) {
    int Nn = Nel + 1;
    std::vector<double> x(Nn);
    double h = L / Nel;
    for (int i = 0; i < Nn; ++i) x[i] = i * h;

    MatrixXd K = MatrixXd::Zero(Nn, Nn);
    VectorXd F = VectorXd::Zero(Nn);

    for (int e = 0; e < Nel; ++e) {
        int i = e, j = e + 1;
        MatrixXd Ke(2,2);
        Ke <<  1.0, -1.0,
              -1.0,  1.0;
        Ke *= (1.0 / h);

        double xm = 0.5 * (x[i] + x[j]);
        double fe = f(xm);
        Eigen::Vector2d Fe;
        Fe << 1.0, 1.0;
        Fe *= (fe * (h / 2.0));

        K.block(i, i, 2, 2) += Ke;
        F.segment(i, 2) += Fe;
    }

    // Dirichlet elimination
    std::vector<int> free;
    free.reserve(Nn - 2);
    for (int i = 1; i < Nn - 1; ++i) free.push_back(i);

    // modify RHS
    F -= K.col(0) * u0;
    F -= K.col(Nn - 1) * uL;

    MatrixXd Kff(free.size(), free.size());
    VectorXd Ff(free.size());
    for (int a = 0; a < (int)free.size(); ++a) {
        Ff(a) = F(free[a]);
        for (int b = 0; b < (int)free.size(); ++b) {
            Kff(a,b) = K(free[a], free[b]);
        }
    }

    VectorXd uf = Kff.fullPivLu().solve(Ff);

    VectorXd u = VectorXd::Zero(Nn);
    u(0) = u0;
    u(Nn - 1) = uL;
    for (int a = 0; a < (int)free.size(); ++a) u(free[a]) = uf(a);

    return {x, u};
}

double f_rhs(double x) {
    // Choose f so exact solution is u=sin(pi x) on [0,1] with u(0)=u(1)=0:
    // -u'' = pi^2 sin(pi x)
    return (M_PI * M_PI) * std::sin(M_PI * x);
}

int main() {
    // ---- Heat FD (theta method) ----
    Grid1D grid(1.0, 101);
    double alpha = 0.05;

    VectorXd u0(grid.N);
    for (int i = 0; i < grid.N; ++i) u0(i) = std::sin(M_PI * grid.x[i]);
    u0(0) = 0.0; u0(grid.N - 1) = 0.0;

    double dt = 0.4 * (grid.dx * grid.dx) / alpha; // stable for explicit, also fine for theta
    int steps = 200;

    std::vector<double> t;
    std::vector<VectorXd> u_hist;
    heat_fd_theta(u0, alpha, grid, dt, steps, 0.5, t, u_hist); // Crank–Nicolson

    std::cout << "Heat FD done. Final time = " << t.back() << "\n";

    // Save final profile to CSV
    {
        std::ofstream ofs("heat_fd_final.csv");
        ofs << "x,u\n";
        const VectorXd& uF = u_hist.back();
        for (int i = 0; i < grid.N; ++i) {
            ofs << grid.x[i] << "," << uF(i) << "\n";
        }
    }

    // ---- Poisson FEM (P1) ----
    auto [x, u] = fem_poisson_p1(40, 1.0, f_rhs, 0.0, 0.0);
    double err_inf = 0.0;
    for (int i = 0; i < (int)x.size(); ++i) {
        double uex = std::sin(M_PI * x[i]);
        err_inf = std::max(err_inf, std::abs(u(i) - uex));
    }
    std::cout << "FEM Poisson max error = " << err_inf << "\n";

    // Save FEM solution
    {
        std::ofstream ofs("poisson_fem.csv");
        ofs << "x,u,uex\n";
        for (int i = 0; i < (int)x.size(); ++i) {
            double uex = std::sin(M_PI * x[i]);
            ofs << x[i] << "," << u(i) << "," << uex << "\n";
        }
    }

    std::cout << "Wrote heat_fd_final.csv and poisson_fem.csv\n";
    return 0;
}
      

File: Chapter19_Lesson3.java

/*
Chapter19_Lesson3.java
Finite Difference (FD) and Finite Element (FE) Approximations (Conceptual Level)

This Java example shows:
  1) 1D heat equation u_t = alpha u_xx via FD in space and theta-method in time.
  2) 1D Poisson -u'' = f(x) via linear FEM (Galerkin) on uniform mesh.

Recommended linear algebra library:
  - EJML (Efficient Java Matrix Library)
    https://ejml.org/

This code uses EJML's SimpleMatrix for clarity (dense matrices). For large N, use EJML sparse.

Compile/run (example with Maven/Gradle):
  Ensure EJML is on the classpath.

Author: (course material generator)
*/

import org.ejml.simple.SimpleMatrix;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Locale;

public class Chapter19_Lesson3 {

    static class Grid1D {
        double L;
        int N;
        double[] x;
        double dx;

        Grid1D(double L, int N) {
            if (N < 3) throw new IllegalArgumentException("N must be >= 3");
            this.L = L; this.N = N;
            this.x = new double[N];
            for (int i = 0; i < N; i++) x[i] = (L * i) / (N - 1);
            this.dx = x[1] - x[0];
        }
    }

    // Dense Laplacian on interior nodes with Dirichlet boundaries.
    static SimpleMatrix laplacianDirichlet1D(int N, double dx) {
        int n = N - 2;
        SimpleMatrix A = new SimpleMatrix(n, n);
        double s = 1.0 / (dx * dx);
        for (int i = 0; i < n; i++) {
            A.set(i, i, -2.0 * s);
            if (i - 1 >= 0) A.set(i, i - 1, 1.0 * s);
            if (i + 1 < n)  A.set(i, i + 1, 1.0 * s);
        }
        return A;
    }

    static double[] heatFDTheta(
            double[] u0Full,
            double alpha,
            Grid1D grid,
            double dt,
            int steps,
            double theta
    ) {
        int N = grid.N;
        int n = N - 2;

        SimpleMatrix A = laplacianDirichlet1D(N, grid.dx);
        SimpleMatrix I = SimpleMatrix.identity(n);

        SimpleMatrix LHS = I.minus(A.scale(theta * dt * alpha));
        SimpleMatrix RHS = I.plus(A.scale((1.0 - theta) * dt * alpha));

        double uL = u0Full[0];
        double uR = u0Full[N - 1];

        // boundary vector bc for interior equation
        double s = 1.0 / (grid.dx * grid.dx);
        SimpleMatrix bc = new SimpleMatrix(n, 1);
        bc.set(0, 0, uL * s);
        bc.set(n - 1, 0, uR * s);

        SimpleMatrix v = new SimpleMatrix(n, 1);
        for (int i = 0; i < n; i++) v.set(i, 0, u0Full[i + 1]);

        for (int k = 0; k < steps; k++) {
            SimpleMatrix rhs = RHS.mult(v).plus(bc.scale(dt * alpha));
            v = LHS.solve(rhs); // dense solve
        }

        double[] uFinal = new double[N];
        uFinal[0] = uL;
        uFinal[N - 1] = uR;
        for (int i = 0; i < n; i++) uFinal[i + 1] = v.get(i, 0);
        return uFinal;
    }

    static class FemResult {
        double[] x;
        double[] u;
        FemResult(double[] x, double[] u) { this.x = x; this.u = u; }
    }

    // Linear FEM for -u'' = f(x) on (0,L), Dirichlet u(0)=u0, u(L)=uL.
    static FemResult femPoissonP1(int Nel, double L, java.util.function.DoubleUnaryOperator f, double u0, double uL) {
        int Nn = Nel + 1;
        double[] x = new double[Nn];
        double h = L / Nel;
        for (int i = 0; i < Nn; i++) x[i] = i * h;

        SimpleMatrix K = new SimpleMatrix(Nn, Nn);
        SimpleMatrix F = new SimpleMatrix(Nn, 1);

        for (int e = 0; e < Nel; e++) {
            int i = e, j = e + 1;
            SimpleMatrix Ke = new SimpleMatrix(2, 2);
            Ke.set(0,0,  1.0/h); Ke.set(0,1, -1.0/h);
            Ke.set(1,0, -1.0/h); Ke.set(1,1,  1.0/h);

            double xm = 0.5 * (x[i] + x[j]);
            double fe = f.applyAsDouble(xm);
            double load = fe * (h / 2.0);

            // assemble
            K.set(i, i, K.get(i,i) + Ke.get(0,0));
            K.set(i, j, K.get(i,j) + Ke.get(0,1));
            K.set(j, i, K.get(j,i) + Ke.get(1,0));
            K.set(j, j, K.get(j,j) + Ke.get(1,1));

            F.set(i, 0, F.get(i,0) + load);
            F.set(j, 0, F.get(j,0) + load);
        }

        // Dirichlet elimination: solve for interior nodes 1..Nn-2
        int n = Nn - 2;
        SimpleMatrix Kff = new SimpleMatrix(n, n);
        SimpleMatrix Ff  = new SimpleMatrix(n, 1);

        for (int a = 0; a < n; a++) {
            int I = a + 1;
            double rhs = F.get(I, 0) - K.get(I, 0) * u0 - K.get(I, Nn - 1) * uL;
            Ff.set(a, 0, rhs);
            for (int b = 0; b < n; b++) {
                int J = b + 1;
                Kff.set(a, b, K.get(I, J));
            }
        }

        SimpleMatrix uf = Kff.solve(Ff);

        double[] u = new double[Nn];
        u[0] = u0;
        u[Nn - 1] = uL;
        for (int a = 0; a < n; a++) u[a + 1] = uf.get(a, 0);

        return new FemResult(x, u);
    }

    static double fRhs(double x) {
        // -u'' = pi^2 sin(pi x) gives exact u = sin(pi x) on [0,1]
        return Math.PI * Math.PI * Math.sin(Math.PI * x);
    }

    static void writeCSV(String filename, double[] x, double[] u, double[] uex) throws IOException {
        try (FileWriter fw = new FileWriter(filename)) {
            fw.write("x,u,uex\n");
            for (int i = 0; i < x.length; i++) {
                fw.write(String.format(Locale.US, "%.10f,%.10f,%.10f\n", x[i], u[i], uex[i]));
            }
        }
    }

    public static void main(String[] args) throws IOException {
        Grid1D grid = new Grid1D(1.0, 101);
        double alpha = 0.05;

        double[] u0 = new double[grid.N];
        for (int i = 0; i < grid.N; i++) u0[i] = Math.sin(Math.PI * grid.x[i]);
        u0[0] = 0.0; u0[grid.N - 1] = 0.0;

        double dt = 0.4 * (grid.dx * grid.dx) / alpha;
        int steps = 200;
        double[] uHeatFinal = heatFDTheta(u0, alpha, grid, dt, steps, 0.5);

        // write heat solution
        double[] uexHeat = new double[grid.N];
        for (int i = 0; i < grid.N; i++) uexHeat[i] = uHeatFinal[i]; // no analytic in this demo
        writeCSV("heat_fd_final_java.csv", grid.x, uHeatFinal, uexHeat);

        // FEM Poisson
        FemResult fr = femPoissonP1(40, 1.0, Chapter19_Lesson3::fRhs, 0.0, 0.0);
        double[] uex = new double[fr.x.length];
        double errInf = 0.0;
        for (int i = 0; i < fr.x.length; i++) {
            uex[i] = Math.sin(Math.PI * fr.x[i]);
            errInf = Math.max(errInf, Math.abs(fr.u[i] - uex[i]));
        }
        System.out.println("FEM Poisson max error = " + errInf);
        writeCSV("poisson_fem_java.csv", fr.x, fr.u, uex);

        System.out.println("Wrote heat_fd_final_java.csv and poisson_fem_java.csv");
    }
}
      

File: Chapter19_Lesson3.m

% Chapter19_Lesson3.m
% Finite Difference (FD) and Finite Element (FE) Approximations (Conceptual Level)
%
% This script provides:
%   1) 1D heat equation via FD (method of lines) + theta method in time.
%   2) 1D Poisson -u''=f(x) via linear FEM (Galerkin) on uniform mesh.
%
% MATLAB toolboxes:
%   - Base MATLAB is sufficient for this script.
%   - (Optional) PDE Toolbox offers higher-level FEM routines, but we implement from scratch here.

clear; clc;

%% 1) Heat equation u_t = alpha u_xx on [0,1], Dirichlet u(0)=u(1)=0
L = 1.0;
N = 101;                 % nodes
x = linspace(0, L, N)';
dx = x(2) - x(1);
alpha = 0.05;

u0 = sin(pi*x);
u0(1) = 0.0; u0(end) = 0.0;

% interior size
n = N - 2;

% Laplacian on interior nodes (Dirichlet)
e = ones(n,1);
A = spdiags([e -2*e e], [-1 0 1], n, n) / (dx*dx);

theta = 0.5;             % Crank-Nicolson
dt = 0.4 * (dx*dx) / alpha;
steps = 200;

I = speye(n);
LHS = (I - theta*dt*alpha*A);
RHS = (I + (1-theta)*dt*alpha*A);

uL = u0(1);
uR = u0(end);

% boundary vector for interior equation
bc = zeros(n,1);
bc(1) = bc(1) + uL/(dx*dx);
bc(end) = bc(end) + uR/(dx*dx);

v = u0(2:end-1);

for k = 1:steps
    rhs = RHS*v + dt*alpha*bc;
    v = LHS \ rhs;
end

uHeatFinal = zeros(N,1);
uHeatFinal(1) = uL;
uHeatFinal(end) = uR;
uHeatFinal(2:end-1) = v;

fprintf('Heat FD done. Final time = %.6f\n', steps*dt);

%% 2) Poisson -u'' = f(x) on [0,1], u(0)=u(1)=0, exact solution u=sin(pi x)
Nel = 40;
xf = linspace(0,1,Nel+1)';
h = xf(2) - xf(1);
Nn = Nel + 1;

K = zeros(Nn, Nn);
F = zeros(Nn, 1);

f = @(xx) (pi^2) * sin(pi*xx);

for e = 1:Nel
    i = e;
    j = e + 1;

    Ke = (1/h) * [ 1 -1; -1 1 ];
    xm = 0.5*(xf(i)+xf(j));
    fe = f(xm);
    Fe = fe*(h/2)*[1;1];

    K(i:i+1, i:i+1) = K(i:i+1, i:i+1) + Ke;
    F(i:i+1) = F(i:i+1) + Fe;
end

% Dirichlet eliminate nodes 1 and Nn
free = 2:Nn-1;
u0_bc = 0.0; uL_bc = 0.0;

Fmod = F - K(:,1)*u0_bc - K(:,end)*uL_bc;
Kff = K(free, free);
Ff  = Fmod(free);

uf = Kff \ Ff;

uFEM = zeros(Nn,1);
uFEM(1) = u0_bc; uFEM(end) = uL_bc;
uFEM(free) = uf;

uExact = sin(pi*xf);
errInf = norm(uFEM - uExact, inf);
fprintf('FEM Poisson max error = %.6e\n', errInf);

%% (Optional) plots
figure; plot(x, uHeatFinal, 'LineWidth', 1.5); grid on;
xlabel('x'); ylabel('u'); title('Heat equation (FD space + theta time)');

figure; plot(xf, uFEM, 'o-', xf, uExact, '--', 'LineWidth', 1.5); grid on;
xlabel('x'); ylabel('u'); title('Poisson equation (linear FEM)'); legend('FEM','Exact');

%% Note on Simulink (conceptual):
% If you form the MOL ODE for interior nodes: v_dot = alpha*A*v + alpha*bc,
% you can implement it in Simulink as:
%   - Integrator block for state v
%   - Matrix Gain block with alpha*A feeding back
%   - Constant block alpha*bc as input
% This yields a large state-space model representing the distributed PDE.
      

File: Chapter19_Lesson3.nb

(* Chapter19_Lesson3.nb
   Finite Difference (FD) and Finite Element (FE) Approximations (Conceptual Level)

   This file is written as a Mathematica script with notebook-style section markers.
   You can open it in Wolfram Mathematica; the FrontEnd recognizes these markers.

   Contents:
     1) Heat equation by FD (method of lines) + theta method.
     2) Poisson equation by linear FEM (Galerkin) assembly.
*)

(* ::Section:: *)
(*1. Heat equation: FD in space, theta method in time*)

ClearAll["Global`*"];

L = 1.0;
N = 101;
x = Subdivide[0, L, N - 1];
dx = x[[2]] - x[[1]];
alpha = 0.05;

u0 = Sin[Pi x];
u0[[1]] = 0.0; u0[[-1]] = 0.0;

n = N - 2; (* interior size *)

(* Laplacian on interior nodes with Dirichlet boundaries *)
main = ConstantArray[-2/(dx^2), n];
off  = ConstantArray[ 1/(dx^2), n - 1];
A = SparseArray[
   {Band[{1, 1}] -> main, Band[{1, 2}] -> off, Band[{2, 1}] -> off},
   {n, n}
];

theta = 1/2;  (* Crank–Nicolson *)
dt = 0.4*(dx^2)/alpha;
steps = 200;

I = IdentityMatrix[n, SparseArray];
LHS = I - theta*dt*alpha*A;
RHS = I + (1 - theta)*dt*alpha*A;

uL = u0[[1]];
uR = u0[[-1]];

bc = ConstantArray[0.0, n];
bc[[1]] = bc[[1]] + uL/(dx^2);
bc[[-1]] = bc[[-1]] + uR/(dx^2);
bc = SparseArray[bc];

v = u0[[2 ;; -2]];

Do[
  rhs = RHS.v + dt*alpha*bc;
  v = LinearSolve[LHS, rhs];
, {k, 1, steps}];

uHeatFinal = Join[{uL}, v, {uR}];

Print["Heat FD done. Final time = ", steps*dt];

(* Optional: plot *)
ListLinePlot[Transpose[{x, uHeatFinal}], PlotRange -> All, AxesLabel -> {"x", "u"},
 PlotLabel -> "Heat equation (FD space + theta time)"];

(* ::Section:: *)
(*2. Poisson equation: linear FEM (P1) on uniform mesh*)

Nel = 40;
xf = Subdivide[0, 1, Nel];
h = xf[[2]] - xf[[1]];
Nn = Nel + 1;

f[x_] := Pi^2 Sin[Pi x]; (* exact u = Sin[Pi x] *)

K = ConstantArray[0.0, {Nn, Nn}];
F = ConstantArray[0.0, Nn];

Do[
  i = e; j = e + 1;
  Ke = (1/h) { {1, -1}, {-1, 1} };
  xm = (xf[[i]] + xf[[j]])/2;
  Fe = f[xm]*(h/2)*{1, 1};

  K[[{i, j}, {i, j}]] += Ke;
  F[[{i, j}]] += Fe;
, {e, 1, Nel}];

(* Dirichlet BC u(0)=u(1)=0, eliminate first and last nodes *)
free = Range[2, Nn - 1];

u0bc = 0.0; uLbc = 0.0;
Fmod = F - K[[All, 1]]*u0bc - K[[All, -1]]*uLbc;

Kff = K[[free, free]];
Ff  = Fmod[[free]];

uf = LinearSolve[Kff, Ff];

uFEM = ConstantArray[0.0, Nn];
uFEM[[1]] = u0bc; uFEM[[-1]] = uLbc;
uFEM[[free]] = uf;

uExact = Sin[Pi xf];
errInf = Max[Abs[uFEM - uExact]];
Print["FEM Poisson max error = ", errInf];

ListLinePlot[{Transpose[{xf, uFEM}], Transpose[{xf, uExact}]},
 PlotLegends -> {"FEM", "Exact"}, PlotRange -> All, AxesLabel -> {"x", "u"},
 PlotLabel -> "Poisson equation (linear FEM)"];
      

Problem 1 (Second-derivative stencil and truncation error): Starting from Taylor expansions of \( u(x_i \pm h) \), derive \( u_{xx}(x_i) \approx (u_{i+1}-2u_i+u_{i-1})/h^2 \) and show truncation error \( \mathcal{O}(h^2) \).

Solution: Add the Taylor series for \( u(x_i+h) \) and \( u(x_i-h) \):

\[ u(x_i+h)+u(x_i-h) = 2u(x_i)+h^2u_{xx}(x_i)+\frac{h^4}{12}u_{xxxx}(x_i)+\mathcal{O}(h^6). \]

Rearrange to obtain:

\[ \frac{u(x_i+h)-2u(x_i)+u(x_i-h)}{h^2} = u_{xx}(x_i)+\frac{h^2}{12}u_{xxxx}(x_i)+\mathcal{O}(h^4). \]

Problem 2 (FTCS stability for heat): For FTCS: \( u_i^{k+1}=u_i^k+r(u_{i+1}^k-2u_i^k+u_{i-1}^k) \), \( r=\alpha\Delta t/h^2 \). Using a Fourier mode \( u_i^k=g^k e^{\mathrm{i}\xi i} \), derive the amplification factor and show \( r \le 1/2 \) suffices for stability.

Solution: Substitute the mode:

\[ g = 1 + r\left(e^{\mathrm{i}\xi}-2+e^{-\mathrm{i}\xi}\right) = 1 - 4r\sin^2\left(\frac{\xi}{2}\right). \]

Since \( \sin^2(\xi/2)\in[0,1] \), the minimum of \( g \) is \( 1-4r \). Enforcing \( |g|\le 1 \) for all modes implies \( 1-4r \ge -1 \), hence \( r \le 1/2 \).

Problem 3 (Local FEM stiffness): On an element \( [x_e,x_{e+1}] \) of length \( h_e \), P1 basis derivatives are \( \varphi_1'=-1/h_e \), \( \varphi_2'=1/h_e \). Compute \( (K_e)_{ij}=\int_{x_e}^{x_{e+1} }\varphi_i'\varphi_j'\,dx \).

Solution: Derivatives are constant, so

\[ (K_e)_{ij} = (\varphi_i'\varphi_j')\int_{x_e}^{x_{e+1} }1\,dx = (\varphi_i'\varphi_j')h_e. \]

Substituting yields \( (K_e)_{11}=1/h_e \), \( (K_e)_{22}=1/h_e \), \( (K_e)_{12}=(K_e)_{21}=-1/h_e \).

Problem 4 (FEM semi-discrete heat system): Starting from \( \int_0^L u_t v\,dx + \alpha\int_0^L u'v'\,dx = 0 \), insert \( u_h(x,t)=\sum_i U_i(t)\varphi_i(x) \) and test with \( v=\varphi_j \). Derive \( \mathbf{M}\dot{\mathbf{U} }+\alpha\mathbf{K}\mathbf{U}=\mathbf{0} \).

Solution:

\[ \sum_i \dot{U}_i(t)\int_0^L \varphi_i\varphi_j\,dx + \alpha\sum_i U_i(t)\int_0^L \varphi_i'\varphi_j'\,dx = 0. \]

Define \( M_{ji}=\int_0^L \varphi_i\varphi_j\,dx \) and \( K_{ji}=\int_0^L \varphi_i'\varphi_j'\,dx \), giving the matrix form.

Problem 5 (FD vs FEM for 1D Poisson on uniform mesh): On a uniform mesh, show the FEM stiffness has the same tridiagonal pattern as the FD Laplacian (up to scaling). Explain the conceptual difference in how the right-hand side is formed.

Solution: Assembly of \( (1/h)\begin{bmatrix}1&-1\\-1&1\end{bmatrix} \) yields a global tridiagonal matrix with diagonal \( 2/h \) and off-diagonals \( -1/h \) (Dirichlet-adjusted). FD yields diagonal \( 2/h^2 \), off-diagonals \( -1/h^2 \). The pattern matches, but FEM forms \( F_j=\int_0^L f\varphi_j\,dx \) (quadrature/integral), whereas FD often samples \( f \) pointwise.

Spatial discretization converts PDEs into high-dimensional state-space dynamics. FD constructs stencil operators with quantifiable truncation error and stability restrictions (often \( \Delta t=\mathcal{O}(h^2) \) for diffusion when explicit). FEM derives from weak forms and yields sparse matrix ODEs involving mass and stiffness matrices, with quasi-optimal error bounds in energy norms. Both paradigms supply the computational bridge from distributed physics to simulation and control-oriented models.