In earlier chapters, most models were lumped-parameter: the system state was described by finitely many variables (e.g., displacement of a mass, voltage of a capacitor, temperature of a lump). In a distributed-parameter system, the state depends on space and time, so the state is a field.

\( u(x,t) \): a spatially distributed state (temperature, concentration, displacement, etc.) over a domain and time.

Typical one-dimensional domains are rods, strings, pipelines, and channels, with \( 0 < x < L \) and \( t > 0 \). The canonical PDEs in this lesson are:

\[ \text{Heat / Diffusion (parabolic):}\quad \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \]

\[ \text{Wave (hyperbolic):}\quad \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \]

Here \( \alpha \) is a diffusivity (thermal or mass diffusivity) and \( c \) is a wave propagation speed. These PDEs are the distributed analogs of first- and second-order ODE dynamics from lumped systems.

flowchart TD
  A["Physical domain (rod / string / channel)"] --> B["Choose field variable u(x,t)"]
  B --> C["Conservation law or force balance on small element"]
  C --> D["Constitutive law (Fourier / Fick / tension law)"]
  D --> E["Take limit as element size shrinks"]
  E --> F["PDE + initial conditions + boundary conditions"]
  F --> G["Analyze: energy, stability, uniqueness"]
  G --> H["Simulate or discretize for computation"]
        

The heat equation and diffusion equation share the same mathematical structure. Only the physical meaning of the state and parameters changes.

2.1 Diffusion Equation from Mass Conservation + Fick's Law

Let \( C(x,t) \) be concentration (mass per unit volume) in a 1D medium with constant cross-sectional area \( A \). Consider a small interval \( [x, x+\Delta x] \). Mass conservation gives:

\[ A \Delta x \frac{\partial C}{\partial t} = A J(x,t) - A J(x+\Delta x,t) \]

where \( J(x,t) \) is the diffusive flux (mass flow per unit area). Using Fick's law:

\[ J(x,t) = -D \frac{\partial C}{\partial x} \]

substitute into conservation:

\[ \Delta x \frac{\partial C}{\partial t} = D \left( \frac{\partial C}{\partial x}(x+\Delta x,t) - \frac{\partial C}{\partial x}(x,t)\right) \]

Divide by \( \Delta x \) and take the limit \( \Delta x \to 0 \) (written in text intentionally; no arrow notation inside MathJax):

\[ \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} \]

This is the one-dimensional diffusion equation.

2.2 Heat Equation from Energy Conservation + Fourier's Law

Let \( T(x,t) \) be temperature in a rod with density \( \rho \), specific heat \( c_p \), conductivity \( k \), and cross-sectional area \( A \). Thermal energy in \( [x, x+\Delta x] \) is:

\[ \rho c_p A \Delta x \, T(x,t) \]

Energy conservation (with no internal heat source yet) gives:

\[ \rho c_p A \Delta x \frac{\partial T}{\partial t} = A q(x,t) - A q(x+\Delta x,t) \]

where \( q(x,t) \) is heat flux. Fourier's law:

\[ q(x,t) = -k \frac{\partial T}{\partial x} \]

Substituting and simplifying:

\[ \rho c_p \frac{\partial T}{\partial t} = k \frac{\partial^2 T}{\partial x^2} \]

Define thermal diffusivity \( \alpha = \dfrac{k}{\rho c_p} \). Then:

\[ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} \]

2.3 Source Terms

If a distributed source exists (heat generation or species generation), the PDE becomes:

\[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} + s(x,t) \]

In 3D notation, the same physics is expressed with the Laplacian:

\[ \frac{\partial u}{\partial t} = \alpha \nabla^2 u + s(\mathbf{x},t) \]

Consider a taut string with uniform linear mass density \( \mu \) and constant tension magnitude \( \mathcal{T} \). Let \( y(x,t) \) denote transverse displacement. For small slopes, the vertical force on a small string segment is approximated by the difference of tension components at the ends.

Newton's second law on \( [x, x+\Delta x] \) gives:

\[ \mu \Delta x \frac{\partial^2 y}{\partial t^2} = \mathcal{T}\sin\theta(x+\Delta x,t) - \mathcal{T}\sin\theta(x,t) \]

For small slopes, \( \sin\theta \approx \tan\theta = \dfrac{\partial y}{\partial x} \), so:

\[ \mu \Delta x \frac{\partial^2 y}{\partial t^2} = \mathcal{T}\left( \frac{\partial y}{\partial x}(x+\Delta x,t) - \frac{\partial y}{\partial x}(x,t)\right) \]

Divide by \( \Delta x \) and take the limit:

\[ \mu \frac{\partial^2 y}{\partial t^2} = \mathcal{T}\frac{\partial^2 y}{\partial x^2} \]

Defining the wave speed \( c = \sqrt{\mathcal{T}/\mu} \), we obtain:

\[ \frac{\partial^2 y}{\partial t^2} = c^2 \frac{\partial^2 y}{\partial x^2} \]

Unlike the heat/diffusion equation, the wave equation is second order in time and stores kinetic and potential energy, which is why it behaves like a distributed version of undamped oscillatory dynamics.

A PDE model is incomplete without initial and boundary conditions. The type of PDE determines how many initial conditions are needed in time.

4.1 Initial Conditions

For the heat/diffusion equation (first order in time), one initial condition is required:

\[ u(x,0) = u_0(x) \]

For the wave equation (second order in time), two initial conditions are required:

\[ u(x,0) = u_0(x), \qquad \frac{\partial u}{\partial t}(x,0) = v_0(x) \]

4.2 Boundary Conditions

At each spatial boundary (e.g., \( x=0 \) and \( x=L \)), common conditions are:

Dirichlet (prescribed state):

\[ u(0,t) = g_0(t), \qquad u(L,t) = g_L(t) \]

Neumann (prescribed flux / gradient):

\[ \frac{\partial u}{\partial x}(0,t) = h_0(t), \qquad \frac{\partial u}{\partial x}(L,t) = h_L(t) \]

Robin (convective / mixed):

\[ -k\frac{\partial T}{\partial x}(0,t) = h\big(T(0,t)-T_{\infty}\big) \]

Robin conditions are especially important in heat transfer, where the boundary exchanges heat with an ambient environment.

4.3 Physical Interpretation in System Dynamics

Boundary conditions play the same role as interconnection constraints in lumped models. A clamped string, insulated rod, and fixed-temperature rod differ primarily in their boundary conditions, even if the PDE in the interior is identical.

Energy methods are central in System Dynamics because they connect the differential equations to physical storage and dissipation. The same philosophy extends to PDEs.

5.1 Heat/Diffusion Equation Energy Decay

Consider the homogeneous heat equation on \( 0 < x < L \):

\[ u_t = \alpha u_{xx}, \qquad \alpha > 0 \]

Define the quadratic energy-like functional:

\[ E_h(t) = \frac{1}{2}\int_0^L u(x,t)^2 \, dx \]

Differentiate with respect to time:

\[ \frac{dE_h}{dt} = \int_0^L u \, u_t \, dx = \alpha \int_0^L u \, u_{xx} \, dx \]

Integration by parts yields:

\[ \frac{dE_h}{dt} = \alpha \left[u u_x\right]_0^L - \alpha \int_0^L u_x^2 \, dx \]

Under homogeneous Dirichlet conditions \( u(0,t)=u(L,t)=0 \), or homogeneous Neumann conditions with compatible boundary term, the boundary term vanishes. Therefore:

\[ \frac{dE_h}{dt} = - \alpha \int_0^L u_x^2 \, dx \le 0 \]

This proves monotonic decay of the energy-like norm. Physically, diffusion smooths gradients and dissipates spatial variation.

5.2 Uniqueness for the Heat Equation (Energy Method)

Suppose \( u_1 \) and \( u_2 \) satisfy the same heat equation, same initial condition, and same homogeneous Dirichlet boundary conditions. Let \( w = u_1 - u_2 \). Then:

\[ w_t = \alpha w_{xx}, \qquad w(x,0)=0, \qquad w(0,t)=w(L,t)=0 \]

Applying the previous energy identity to \( w \):

\[ \frac{d}{dt}\left(\frac{1}{2}\int_0^L w^2 dx\right) = -\alpha \int_0^L w_x^2 dx \le 0 \]

Since the initial energy is zero, it remains zero for all time:

\[ \int_0^L w(x,t)^2 dx = 0 \]

Hence \( w(x,t)=0 \) almost everywhere, so \( u_1=u_2 \). This establishes uniqueness.

5.3 Wave Equation Energy Conservation

For the homogeneous wave equation:

\[ u_{tt} = c^2 u_{xx} \]

define the wave energy:

\[ E_w(t) = \frac{1}{2}\int_0^L \left(u_t^2 + c^2 u_x^2\right)dx \]

Differentiate:

\[ \frac{dE_w}{dt} = \int_0^L u_t u_{tt}\,dx + c^2 \int_0^L u_x u_{xt}\,dx \]

Substitute \( u_{tt}=c^2 u_{xx} \):

\[ \frac{dE_w}{dt} = c^2\int_0^L u_t u_{xx}\,dx + c^2\int_0^L u_x u_{xt}\,dx \]

Integrate the first integral by parts in space:

\[ \int_0^L u_t u_{xx}\,dx = \left[u_t u_x\right]_0^L - \int_0^L u_{xt}u_x\,dx \]

The two volume integrals cancel, so:

\[ \frac{dE_w}{dt} = c^2\left[u_t u_x\right]_0^L \]

Under fixed-end (Dirichlet) or insulated-type compatible boundary conditions that make the boundary power term zero:

\[ \frac{dE_w}{dt} = 0 \]

Therefore wave dynamics conserves energy (in the ideal lossless case), in contrast to heat/diffusion dynamics.

In state-space modeling, the state is usually a finite vector \( \mathbf{x}(t) \). For PDEs, the state is a function over space. A useful mental model is:

\( x(t) \): not a finite vector, but a function \( u(\cdot,t) \) in a function space (for example, square-integrable functions over the domain).

This is why distributed-parameter systems are often called infinite-dimensional systems. The PDE defines the time evolution operator on this function-valued state.

6.1 Semi-Discrete Preview (without full numerical theory yet)

If we discretize space into grid points (formal details are developed in Lesson 3), the PDE becomes a large ODE system. For the heat equation:

\[ \dot{\mathbf{u}}(t) = \alpha \mathbf{D}_{xx}\mathbf{u}(t) + \mathbf{s}(t) \]

and for the wave equation:

\[ \ddot{\mathbf{u}}(t) = c^2 \mathbf{D}_{xx}\mathbf{u}(t) \]

where \( \mathbf{D}_{xx} \) is a matrix approximation of the second spatial derivative. This explains why many simulation tools (including MATLAB/Simulink) can model PDEs after spatial discretization.

flowchart TD
  A["PDE state u(x,t)"] --> B["Spatial grid x1...xN"]
  B --> C["Approximate second derivative by matrix Dxx"]
  C --> D["Heat: du/dt = alpha Dxx u + s"]
  C --> E["Wave: d2u/dt2 = c^2 Dxx u"]
  D --> F["Use ODE solvers / state-space tools"]
  E --> F
        

6.2 Nondimensionalization and Time Scales

Nondimensionalization is important in system dynamics because it reveals dominant effects and characteristic scales. For heat/diffusion in a rod of length \( L \), define:

\[ \xi = \frac{x}{L}, \qquad \theta = \frac{\alpha t}{L^2} \]

Then the diffusion equation becomes:

\[ \frac{\partial u}{\partial \theta} = \frac{\partial^2 u}{\partial \xi^2} \]

showing the characteristic diffusion time scale:

\[ t_d = \frac{L^2}{\alpha} \]

For wave propagation, a characteristic time scale is:

\[ t_w = \frac{L}{c} \]

These two scales are fundamentally different: diffusion time scales with \( L^2 \), while wave time scales with \( L \).

The following script implements explicit finite-difference schemes for the 1D heat/diffusion and wave equations. We use:

\( r = \alpha \Delta t / \Delta x^2 \): diffusion stability ratio (FTCS requires \( r \le 1/2 \)).

\( s = c \Delta t / \Delta x \): wave CFL number (central scheme typically requires \( s \le 1 \)).

File: Chapter19_Lesson1.py

# Chapter19_Lesson1.py
# 1D heat/diffusion and wave equations (explicit finite differences)

import numpy as np
import matplotlib.pyplot as plt

# --- Heat / diffusion: u_t = alpha * u_xx ---
L = 1.0
alpha = 0.2
Nx = 81
x = np.linspace(0.0, L, Nx)
dx = x[1] - x[0]
r = 0.45  # stability: r <= 0.5
dt = r * dx * dx / alpha
T = 0.25
Nt = int(T / dt)

u = np.sin(np.pi * x) + 0.2 * np.sin(3 * np.pi * x)
u[0] = 0.0
u[-1] = 0.0
heat_snaps = [u.copy()]
heat_times = [0.0]

for n in range(Nt):
    un = u.copy()
    u[1:-1] = un[1:-1] + alpha * dt / dx**2 * (un[2:] - 2 * un[1:-1] + un[:-2])
    u[0] = 0.0
    u[-1] = 0.0
    if n % max(1, Nt // 4) == 0:
        heat_snaps.append(u.copy())
        heat_times.append((n + 1) * dt)

# --- Wave: u_tt = c^2 * u_xx ---
c = 1.0
Nxw = 201
xw = np.linspace(0.0, L, Nxw)
dxw = xw[1] - xw[0]
s = 0.95  # CFL: s <= 1
dtw = s * dxw / c
Tw = 1.0
Ntw = int(Tw / dtw)

u0 = np.exp(-180 * (xw - 0.35) ** 2)  # initial displacement
v0 = np.zeros_like(xw)                 # initial velocity
u0[0] = 0.0
u0[-1] = 0.0

u_prev = u0.copy()
u_curr = u0.copy()
u_curr[1:-1] = (
    u0[1:-1]
    + dtw * v0[1:-1]
    + 0.5 * (c * dtw / dxw) ** 2 * (u0[2:] - 2 * u0[1:-1] + u0[:-2])
)
u_curr[0] = 0.0
u_curr[-1] = 0.0

wave_snaps = [u_prev.copy(), u_curr.copy()]
wave_times = [0.0, dtw]

for n in range(1, Ntw):
    u_next = np.zeros_like(u_curr)
    u_next[1:-1] = (
        2 * u_curr[1:-1] - u_prev[1:-1]
        + (c * dtw / dxw) ** 2 * (u_curr[2:] - 2 * u_curr[1:-1] + u_curr[:-2])
    )
    u_prev, u_curr = u_curr, u_next
    if n % max(1, Ntw // 4) == 0:
        wave_snaps.append(u_curr.copy())
        wave_times.append((n + 1) * dtw)

# Plot snapshots
plt.figure(figsize=(8, 4))
for k, t in enumerate(heat_times):
    plt.plot(x, heat_snaps[k], label=f"t={t:.3f}")
plt.title("Heat / Diffusion Equation")
plt.xlabel("x")
plt.ylabel("u(x,t)")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()

plt.figure(figsize=(8, 4))
for k, t in enumerate(wave_times[:5]):
    plt.plot(xw, wave_snaps[k], label=f"t={t:.3f}")
plt.title("Wave Equation")
plt.xlabel("x")
plt.ylabel("u(x,t)")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()
      

This script intentionally uses transparent finite-difference loops and array operations so students can see the PDE-to-discrete-state connection directly before formal numerical-analysis coverage in later lessons.

This C++ implementation writes CSV outputs for post-processing (e.g., plotting in Python, MATLAB, or spreadsheet tools). It is suitable for fast simulation or for embedding into larger control-oriented software pipelines.

File: Chapter19_Lesson1.cpp

// Chapter19_Lesson1.cpp
// 1D heat/diffusion and wave equations (explicit finite differences)
// Compile: g++ -O2 -std=c++17 Chapter19_Lesson1.cpp -o ch19_l1

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

int main() {
    // Heat / diffusion equation: u_t = alpha u_xx
    const double L = 1.0, alpha = 0.2;
    const int Nx = 81;
    const double dx = L / (Nx - 1);
    const double r = 0.45;                 // r <= 0.5
    const double dt = r * dx * dx / alpha;
    const int Nt = static_cast<int>(0.25 / dt);

    std::vector<double> x(Nx), u(Nx), un(Nx);
    for (int i = 0; i < Nx; ++i) {
        x[i] = i * dx;
        u[i] = std::sin(M_PI * x[i]) + 0.2 * std::sin(3.0 * M_PI * x[i]);
    }
    u.front() = 0.0; u.back() = 0.0;

    std::ofstream heatFile("Chapter19_Lesson1_heat.csv");
    heatFile << "step,x,u\n";
    for (int n = 0; n <= Nt; ++n) {
        for (int i = 0; i < Nx; ++i) { heatFile << n << "," << x[i] << "," << u[i] << "\n"; }
        un = u;
        for (int i = 1; i < Nx - 1; ++i) {
            u[i] = un[i] + alpha * dt / (dx * dx) * (un[i + 1] - 2.0 * un[i] + un[i - 1]);
        }
        u.front() = 0.0; u.back() = 0.0;
    }

    // Wave equation: u_tt = c^2 u_xx
    const double c = 1.0;
    const int Nxw = 201;
    const double dxw = L / (Nxw - 1);
    const double s = 0.95;                 // s <= 1
    const double dtw = s * dxw / c;
    const int Ntw = static_cast<int>(1.0 / dtw);

    std::vector<double> xw(Nxw), up(Nxw), uc(Nxw), unext(Nxw, 0.0);
    for (int i = 0; i < Nxw; ++i) {
        xw[i] = i * dxw;
        up[i] = std::exp(-180.0 * (xw[i] - 0.35) * (xw[i] - 0.35));
    }
    up.front() = 0.0; up.back() = 0.0;

    uc = up;
    double cfl2 = (c * dtw / dxw) * (c * dtw / dxw);
    for (int i = 1; i < Nxw - 1; ++i) {
        uc[i] = up[i] + 0.5 * cfl2 * (up[i + 1] - 2.0 * up[i] + up[i - 1]);
    }
    uc.front() = 0.0; uc.back() = 0.0;

    std::ofstream waveFile("Chapter19_Lesson1_wave.csv");
    waveFile << "step,x,u\n";
    for (int n = 0; n <= Ntw; ++n) {
        for (int i = 0; i < Nxw; ++i) { waveFile << n << "," << xw[i] << "," << uc[i] << "\n"; }
        for (int i = 1; i < Nxw - 1; ++i) {
            unext[i] = 2.0 * uc[i] - up[i] + cfl2 * (uc[i + 1] - 2.0 * uc[i] + uc[i - 1]);
        }
        unext.front() = 0.0; unext.back() = 0.0;
        up = uc;
        uc = unext;
    }

    std::cout << "Saved Chapter19_Lesson1_heat.csv and Chapter19_Lesson1_wave.csv\n";
    return 0;
}
      

The C++ version is intentionally compact and uses standard-library containers only, which is useful for teaching the data layout of PDE state vectors and for later integration with control or estimation modules.

Java is often used for engineering backends, digital twins, and teaching tools. The following code mirrors the same finite-difference logic and writes CSV outputs.

File: Chapter19_Lesson1.java

// Chapter19_Lesson1.java
// 1D heat/diffusion and wave equations (explicit finite differences)

import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;

public class Chapter19_Lesson1 {
    public static void main(String[] args) throws IOException {
        // Heat / diffusion: u_t = alpha u_xx
        double L = 1.0, alpha = 0.2;
        int Nx = 81;
        double dx = L / (Nx - 1);
        double r = 0.45;                     // r <= 0.5
        double dt = r * dx * dx / alpha;
        int Nt = (int)(0.25 / dt);

        double[] x = new double[Nx];
        double[] u = new double[Nx];
        double[] un = new double[Nx];
        for (int i = 0; i < Nx; i++) {
            x[i] = i * dx;
            u[i] = Math.sin(Math.PI * x[i]) + 0.2 * Math.sin(3.0 * Math.PI * x[i]);
        }
        u[0] = 0.0; u[Nx - 1] = 0.0;

        try (PrintWriter heatFile = new PrintWriter(new FileWriter("Chapter19_Lesson1_heat_java.csv"))) {
            heatFile.println("step,x,u");
            for (int n = 0; n <= Nt; n++) {
                for (int i = 0; i < Nx; i++) {
                    heatFile.println(n + "," + x[i] + "," + u[i]);
                    un[i] = u[i];
                }
                for (int i = 1; i < Nx - 1; i++) {
                    u[i] = un[i] + alpha * dt / (dx * dx) * (un[i + 1] - 2.0 * un[i] + un[i - 1]);
                }
                u[0] = 0.0; u[Nx - 1] = 0.0;
            }
        }

        // Wave: u_tt = c^2 u_xx
        double c = 1.0;
        int Nxw = 201;
        double dxw = L / (Nxw - 1);
        double s = 0.95;                     // s <= 1
        double dtw = s * dxw / c;
        int Ntw = (int)(1.0 / dtw);

        double[] xw = new double[Nxw];
        double[] up = new double[Nxw];
        double[] uc = new double[Nxw];
        double[] unext = new double[Nxw];
        for (int i = 0; i < Nxw; i++) {
            xw[i] = i * dxw;
            up[i] = Math.exp(-180.0 * (xw[i] - 0.35) * (xw[i] - 0.35));
        }
        up[0] = 0.0; up[Nxw - 1] = 0.0;

        System.arraycopy(up, 0, uc, 0, Nxw);
        double cfl2 = (c * dtw / dxw) * (c * dtw / dxw);
        for (int i = 1; i < Nxw - 1; i++) {
            uc[i] = up[i] + 0.5 * cfl2 * (up[i + 1] - 2.0 * up[i] + up[i - 1]);
        }
        uc[0] = 0.0; uc[Nxw - 1] = 0.0;

        try (PrintWriter waveFile = new PrintWriter(new FileWriter("Chapter19_Lesson1_wave_java.csv"))) {
            waveFile.println("step,x,u");
            for (int n = 0; n <= Ntw; n++) {
                for (int i = 0; i < Nxw; i++) {
                    waveFile.println(n + "," + xw[i] + "," + uc[i]);
                }
                for (int i = 1; i < Nxw - 1; i++) {
                    unext[i] = 2.0 * uc[i] - up[i] + cfl2 * (uc[i + 1] - 2.0 * uc[i] + uc[i - 1]);
                }
                unext[0] = 0.0; unext[Nxw - 1] = 0.0;
                double[] temp = up; up = uc; uc = unext; unext = temp;
            }
        }

        System.out.println("Saved CSV outputs for heat and wave equations.");
    }
}
      

MATLAB is widely used in control engineering courses. The script below uses explicit finite differences and plots solution snapshots. In Simulink, the same semi-discrete equations can be implemented using Integrator blocks after writing the spatially discretized state equations.

File: Chapter19_Lesson1.m

% Chapter19_Lesson1.m
% 1D heat/diffusion and wave equations (explicit finite differences)

clear; clc; close all;

%% Heat / diffusion equation: u_t = alpha u_xx
L = 1.0; alpha = 0.2;
Nx = 81; x = linspace(0, L, Nx); dx = x(2) - x(1);
r = 0.45;                % stability: r <= 0.5
dt = r * dx^2 / alpha;
Nt = floor(0.25 / dt);

u = sin(pi*x) + 0.2*sin(3*pi*x);
u(1) = 0; u(end) = 0;

figure; hold on;
for n = 0:Nt
    if mod(n, max(1, floor(Nt/4))) == 0
        plot(x, u, 'DisplayName', sprintf('t = %.3f', n*dt));
    end
    un = u;
    u(2:end-1) = un(2:end-1) + alpha*dt/dx^2 * (un(3:end) - 2*un(2:end-1) + un(1:end-2));
    u(1) = 0; u(end) = 0;
end
grid on; xlabel('x'); ylabel('u(x,t)'); title('Heat / Diffusion Equation'); legend;

%% Wave equation: u_tt = c^2 u_xx
c = 1.0;
Nxw = 201; xw = linspace(0, L, Nxw); dxw = xw(2) - xw(1);
s = 0.95;                % CFL: s <= 1
dtw = s * dxw / c;
Ntw = floor(1.0 / dtw);

up = exp(-180*(xw - 0.35).^2); up(1) = 0; up(end) = 0;
uc = up;
cfl2 = (c*dtw/dxw)^2;
uc(2:end-1) = up(2:end-1) + 0.5*cfl2 * (up(3:end) - 2*up(2:end-1) + up(1:end-2));

figure; hold on;
for n = 0:Ntw
    if mod(n, max(1, floor(Ntw/4))) == 0
        plot(xw, uc, 'DisplayName', sprintf('t = %.3f', n*dtw));
    end
    unext = zeros(size(uc));
    unext(2:end-1) = 2*uc(2:end-1) - up(2:end-1) + cfl2 * (uc(3:end) - 2*uc(2:end-1) + uc(1:end-2));
    up = uc; uc = unext;
end
grid on; xlabel('x'); ylabel('u(x,t)'); title('Wave Equation'); legend;

%% Simulink note
% The same semi-discrete models can be built in Simulink by stacking Integrator
% blocks for the state-space form obtained after spatial discretization.
      

Mathematica is excellent for PDE prototyping and visualization with NDSolveValue. The following Wolfram Language code solves the same heat/diffusion and wave PDEs.

File: Chapter19_Lesson1.nb

(* Chapter19_Lesson1.nb *)
(* Wolfram Language code for a Mathematica notebook *)

ClearAll["Global`*"];

alpha = 0.2;
heat = NDSolveValue[
  {
    D[u[x, t], t] == alpha D[u[x, t], {x, 2}],
    u[0, t] == 0,
    u[1, t] == 0,
    u[x, 0] == Sin[Pi x] + 0.2 Sin[3 Pi x]
  },
  u,
  {x, 0, 1}, {t, 0, 0.25}
];

Plot[
  Evaluate@Table[heat[x, tt], {tt, {0.00, 0.05, 0.10, 0.20, 0.25}}],
  {x, 0, 1},
  PlotLegends -> {"0.00", "0.05", "0.10", "0.20", "0.25"},
  PlotLabel -> "Heat / Diffusion Equation"
];

c = 1.0;
wave = NDSolveValue[
  {
    D[w[x, t], {t, 2}] == c^2 D[w[x, t], {x, 2}],
    w[0, t] == 0,
    w[1, t] == 0,
    w[x, 0] == Exp[-180 (x - 0.35)^2],
    Derivative[0, 1][w][x, 0] == 0
  },
  w,
  {x, 0, 1}, {t, 0, 1.0}
];

Plot[
  Evaluate@Table[wave[x, tt], {tt, {0.00, 0.20, 0.40, 0.60, 0.80}}],
  {x, 0, 1},
  PlotLegends -> {"0.00", "0.20", "0.40", "0.60", "0.80"},
  PlotLabel -> "Wave Equation"
];
      

Problem 1 (Derivation of Heat Equation): Starting from energy conservation on a rod segment \( [x, x+\Delta x] \) and Fourier's law \( q=-kT_x \), derive \( T_t = \alpha T_{xx} \) and identify \( \alpha \).

Solution:

Energy in the segment is \( \rho c_p A \Delta x \, T \). Conservation gives \( \rho c_p A \Delta x \, T_t = A q(x,t) - A q(x+\Delta x,t) \). Using Fourier's law:

\[ \rho c_p \Delta x \, T_t = k\left(T_x(x+\Delta x,t)-T_x(x,t)\right) \]

Divide by \( \Delta x \) and take the limit:

\[ \rho c_p T_t = k T_{xx} \]

Therefore:

\[ T_t = \alpha T_{xx}, \qquad \alpha = \frac{k}{\rho c_p} \]

Problem 2 (Heat Equation Energy Decay): For \( u_t=\alpha u_{xx} \) with homogeneous Dirichlet boundary conditions, prove that \( E_h(t)=\dfrac{1}{2}\int_0^L u^2 dx \) is non-increasing.

Solution:

Differentiate and substitute the PDE:

\[ \frac{dE_h}{dt} = \int_0^L u u_t dx = \alpha \int_0^L u u_{xx} dx \]

Integration by parts gives:

\[ \frac{dE_h}{dt} = \alpha [u u_x]_0^L - \alpha \int_0^L u_x^2 dx \]

Since \( u(0,t)=u(L,t)=0 \), the boundary term is zero:

\[ \frac{dE_h}{dt} = -\alpha \int_0^L u_x^2 dx \le 0 \]

Hence the energy-like norm decays monotonically.

Problem 3 (Wave Equation Energy Conservation): For \( u_{tt}=c^2u_{xx} \) with boundary conditions that make \( [u_t u_x]_0^L=0 \), prove conservation of \( E_w(t)=\dfrac{1}{2}\int_0^L (u_t^2+c^2u_x^2)\,dx \).

Solution:

Differentiate:

\[ \frac{dE_w}{dt} = \int_0^L u_t u_{tt} dx + c^2 \int_0^L u_x u_{xt} dx \]

Substitute \( u_{tt}=c^2u_{xx} \) and integrate the first term by parts:

\[ \frac{dE_w}{dt} = c^2 [u_t u_x]_0^L - c^2\int_0^L u_{xt}u_x dx + c^2\int_0^L u_xu_{xt} dx \]

The two integrals cancel, leaving:

\[ \frac{dE_w}{dt} = c^2 [u_t u_x]_0^L = 0 \]

Thus \( E_w(t) \) is constant.

Problem 4 (Characteristic Time Scales): A rod of length \( L \) has thermal diffusivity \( \alpha \). A string of the same length has wave speed \( c \). Derive the characteristic diffusion and wave times and compare their dependence on length.

Solution:

Diffusion equation nondimensionalization with \( \xi=x/L \) and \( \theta=\alpha t/L^2 \) yields the natural time scale:

\[ t_d = \frac{L^2}{\alpha} \]

For the wave equation, the propagation time across the domain is:

\[ t_w = \frac{L}{c} \]

Therefore diffusion slows quadratically with domain size, while wave travel time grows linearly with domain size.

Problem 5 (Explicit Scheme Stability Parameters): For the FTCS heat scheme and the centered wave scheme, identify the nondimensional step-size ratios and state the common stability bounds.

Solution:

For heat/diffusion:

\[ r = \frac{\alpha \Delta t}{\Delta x^2} \]

A standard stability condition for the explicit FTCS scheme in 1D is:

\[ r \le \frac{1}{2} \]

For the wave equation:

\[ s = \frac{c\Delta t}{\Delta x} \]

A standard CFL condition for the centered time-space scheme is:

\[ s \le 1 \]

These are precisely the ratios used in the provided Python, C++, Java, and MATLAB scripts.

Problem 6 (Eigenfunction Check for a Fixed-Fixed String/Rod): Verify that \( \phi_n(x)=\sin(n\pi x/L) \) satisfies homogeneous Dirichlet boundary conditions and is an eigenfunction of the second-derivative operator.

Solution:

Boundary conditions:

\[ \phi_n(0)=\sin(0)=0, \qquad \phi_n(L)=\sin(n\pi)=0 \]

Second derivative:

\[ \phi_n''(x) = -\left(\frac{n\pi}{L}\right)^2 \phi_n(x) \]

Thus \( \phi_n \) is an eigenfunction of the Laplacian in 1D with eigenvalue \( -\left(n\pi/L\right)^2 \). This result will be used directly in the next lesson (separation of variables and modal decomposition).

We introduced distributed-parameter modeling through the heat equation, diffusion equation, and wave equation, derived them from first principles (conservation laws and force balance), and formalized initial and boundary conditions as essential parts of a PDE model. We also proved two core energy results: diffusion dynamics dissipates an energy-like norm, while ideal wave dynamics conserves total energy.

These ideas form the mathematical and physical foundation for the next lesson, where we solve PDEs analytically using separation of variables and modal decomposition.