Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section, you will be able to:

Apply the separation of variables technique to reduce the heat equation PDE to two ODEs
Identify the eigenvalue problem that arises from boundary conditions and explain why only discrete eigenvalues are allowed
Construct the Fourier series solution as a superposition of eigenmodes
Calculate Fourier sine coefficients for various initial conditions using orthogonality
Interpret modal decay rates and explain why higher modes decay faster
Connect Fourier analysis to diffusion models in machine learning
Implement the Fourier series solution in Python

The Big Picture: Why Fourier Series?

"Fourier made the remarkable discovery that any periodic function can be expressed as a sum of sines and cosines." — Richard Feynman

In the previous sections, we derived the heat equation and set up the boundary value problem for a rod with fixed temperatures at both ends. Now we face the central question: How do we actually solve it?

The answer lies in one of the most powerful techniques in mathematical physics: Fourier series. The key insight is that certain special solutions—called eigenmodes or normal modes—have a particularly simple form: they maintain their spatial shape while their amplitude decays exponentially.

The Core Strategy

Instead of solving for the temperature u(x, t) directly, we:

Find all special solutions of the form X(x)·T(t) that satisfy the PDE and boundary conditions
Express the initial condition as a sum (Fourier series) of these special solutions
Each mode then evolves independently according to its own decay rate

This approach transforms a complicated PDE problem into simpler ODE problems, plus a Fourier series expansion. It's a beautiful example of how choosing the right basis can dramatically simplify a problem.

Historical Context

The connection between heat diffusion and Fourier series is not coincidental—Fourier invented his series specifically to solve the heat equation!

1807: Fourier's Revolution

Joseph Fourier presented his work on heat conduction to the French Academy. He claimed that any function could be represented as a sum of sines and cosines—a claim so bold that Lagrange and other mathematicians initially rejected it.

The Physical Insight

Fourier realized that heat flow through a rod could be decomposed into independent "modes," each with its own wavelength and decay rate. This physical intuition preceded the rigorous mathematical theory of Fourier series.

Modern Legacy

Fourier analysis now permeates all of science and engineering: from signal processing and image compression to quantum mechanics and—most recently—diffusion models in generative AI.

Separation of Variables

Separation of variables is the systematic approach to finding the eigenmodes. We assume a solution of the form u(x, t) = X(x)·T(t) and show that this leads to a pair of ordinary differential equations.

📐Separation of Variables: Step-by-Step Derivation

Step 1: Start with the Heat Equation

Step 1 of 8

We want to solve the one-dimensional heat equation on a rod of length L with fixed (zero) temperature at both ends:

PDE: ∂u/∂t = α ∂²u/∂x²
Boundary conditions: u(0, t) = 0 and u(L, t) = 0
Initial condition: u(x, 0) = f(x)

Key Equation

fracpartial upartial t = alpha fracpartial^2 upartial x^2

💡 Key Insight

The key idea: we'll guess that u(x,t) can be written as a product of a function of x alone and a function of t alone.

Summary of the Method

The separation of variables technique reveals that the heat equation has solutions that are products of:

Spatial Part: X_n(x)

X_n(x) = sin(nπx/L)

These are eigenfunctions of the Laplacian ∂²/∂x² that automatically satisfy the boundary conditions u(0,t) = u(L,t) = 0.

Temporal Part: T_n(t)

T_n(t) = e^{−α(nπ/L)²t}

Each mode decays exponentially with a rate proportional to n². Higher modes decay faster.

The Eigenvalue Problem

The spatial equation X'' = −λX with boundary conditions X(0) = X(L) = 0 is a Sturm-Liouville eigenvalue problem. This type of problem has deep mathematical significance.

The Eigenvalue Problem

\frac{d^2 X}{dx^2} = -\lambda X, \quad X(0) = 0, \quad X(L) = 0

Solutions exist only for specific values of λ (the eigenvalues):

\lambda_n = \left(\frac{n\pi}{L}\right)^2, \quad n = 1, 2, 3, \ldots

with corresponding eigenfunctions:

X_n(x) = \sin\left(\frac{n\pi x}{L}\right)

Why Discrete Eigenvalues?

The boundary conditions quantize the allowed solutions. Just as a vibrating guitar string can only sustain certain frequencies, the heat equation with fixed boundaries only permits certain "modes" of spatial variation.

Mode n	Eigenvalue λₙ	Spatial Pattern	Physical Meaning
1	π²/L²	sin(πx/L)	Fundamental mode (half sine wave)
2	4π²/L²	sin(2πx/L)	First harmonic (full sine wave)
3	9π²/L²	sin(3πx/L)	Second harmonic (1.5 sine waves)
n	n²π²/L²	sin(nπx/L)	nth harmonic (n half-waves)

Orthogonality

The eigenfunctions sin(nπx/L) are orthogonal over [0, L]:

∫₀^L sin(mπx/L) sin(nπx/L) dx = (L/2)δ_mn

This orthogonality is the key to computing Fourier coefficients—we can "project" any initial condition onto each eigenmode independently.

The Fourier Series Solution

By the principle of superposition (the heat equation is linear), any sum of eigenmodes is also a solution. The general solution is:

The Fourier Series Solution

u(x, t) = \sum_{n=1}^{\infty} b_n \sin\left(\frac{n\pi x}{L}\right) e^{-\alpha (n\pi/L)^2 t}

b_n:Fourier coefficient (determined by initial condition)

sin(nπx/L):Spatial eigenfunction (fixed shape)

e^{−α(nπ/L)²t}:Temporal decay factor (mode-dependent rate)

What This Solution Tells Us

Decomposition: Any initial temperature distribution can be written as a sum of sine waves
Independent evolution: Each mode evolves independently with its own decay rate
Smoothing: Higher modes (sharp features) decay faster, explaining why heat diffusion "smooths" temperature profiles
Long-time behavior: Eventually only the first mode remains, giving u → b₁ sin(πx/L) e^{−α(π/L)²t}

Computing Fourier Coefficients

The coefficients b_n are determined by the initial condition f(x) = u(x, 0). Using the orthogonality of the sine functions:

Fourier Coefficient Formula

b_n = \frac{2}{L} \int_0^L f(x) \sin\left(\frac{n\pi x}{L}\right) dx

This formula "projects" the initial condition f(x) onto the nth eigenmode, measuring how much of that mode is present in the initial temperature distribution.

📊Fourier Sine Coefficients

Number of Terms: 10

Original Function vs Fourier Approximation

Fourier Sine Coefficients b_n

n	b_n	Formula	\|b_n\|/\|b₁\|
1	0.81057	(8/(1²π²))·sin(1π/2)	100.0%
2	0.00000	0 (even n)	0.0%
3	-0.09006	(8/(3²π²))·sin(3π/2)	11.1%
4	0.00000	0 (even n)	0.0%
5	0.03242	(8/(5²π²))·sin(5π/2)	4.0%
6	0.00000	0 (even n)	0.0%
7	-0.01654	(8/(7²π²))·sin(7π/2)	2.0%
8	0.00000	0 (even n)	0.0%
9	0.01001	(8/(9²π²))·sin(9π/2)	1.2%
10	0.00000	0 (even n)	0.0%

The Fourier Sine Coefficients

For the heat equation with Dirichlet boundary conditions u(0,t) = u(L,t) = 0, we expand the initial condition in a sine series:

f(x) = Σ_n=1^∞ b_n sin(nπx/L)

Computing Coefficients

The coefficients are found using the orthogonality of sine functions:

b_n = (2/L) ∫₀^L f(x) sin(nπx/L) dx

Key Insights

• Decay rate: Higher modes (larger n) typically have smaller coefficients
• Symmetry: Symmetric functions (like triangle) have zero even coefficients
• Gibbs phenomenon: Step functions show oscillations near discontinuities
• Convergence: More terms = better approximation (pointwise, except at jumps)

Example: Triangular Initial Condition

Consider a triangular temperature profile peaked at the center:

f(x) = 2x for x ≤ L/2, f(x) = 2(L − x) for x > L/2

Computing the integral:

b_n = \frac{8}{n^2\pi^2} \sin\left(\frac{n\pi}{2}\right)

This gives:

n	sin(nπ/2)	bₙ
1	1	8/π² ≈ 0.811
2	0	0
3	-1	-8/(9π²) ≈ -0.090
4	0	0
5	1	8/(25π²) ≈ 0.032

Pattern in Coefficients

For symmetric initial conditions (like our triangle), the even coefficients are zero. The odd coefficients decay as 1/n², reflecting the smoothness of the triangle (corners cause slow convergence).

The most profound aspect of the Fourier series solution is how each mode decays at its own rate. The decay rate for mode n is:

\text{Decay rate} = \alpha \lambda_n = \alpha \left(\frac{n\pi}{L}\right)^2

🌊Fourier Mode Decomposition

Number of Modes: 5

Thermal Diffusivity α: 0.010

Time t: 0.00

Show modes

Mode Decay Rates

Each mode n decays as e^−λ_n²αt where λ_n = nπ/L

n=1:b₁ = 0.811100.0% remaining

n=3:b₃ = -0.090100.0% remaining

n=5:b₅ = 0.032100.0% remaining

Key Observations

•Higher modes (larger n) decay faster because λ_n² grows as n²
•At t = 0, all modes contribute to the initial condition
•As t → ∞, only the fundamental mode (n=1) survives
•Increasing α speeds up all decay rates proportionally
•This explains why heat diffusion "smooths out" temperature profiles

The Fourier Series Solution

The temperature distribution u(x, t) is expressed as an infinite sum of modes:

u(x, t) = Σ_n=1^∞ b_n · sin(nπx/L) · e^{−α(nπ/L)²t}

Each term sin(nπx/L) is a spatial eigenfunction of the Laplacian, and e^{−α(nπ/L)²t} describes how that mode's amplitude decays over time.

Why Do Higher Modes Decay Faster?

Mathematical Reason

The decay rate is proportional to λ_n = (nπ/L)², which grows as n². Mode 2 decays 4× faster than mode 1, mode 3 decays 9× faster, and so on.

Physical Reason

Higher modes have shorter wavelengths = steeper temperature gradients. Steeper gradients drive faster heat flow (Fourier's law), so sharp features smooth out quickly.

The Smoothing Effect

This differential decay explains why the heat equation is a smoothing operator:

Sharp corners and discontinuities (high-frequency content) disappear quickly
Broad, smooth features (low-frequency content) persist longer
The final state is the smoothest possible: the fundamental mode

Interactive Solution Explorer

Explore how the heat equation solution evolves for different initial conditions. Watch how sharp features smooth out as higher modes decay.

🔥Fourier Series Solution of the Heat Equation

Fourier Modes: 30

Diffusivity α: 0.020

Time t: 0.00

Space-Time Evolution (Heatmap)

Cold

Hot

What You're Seeing

• The dashed line shows the initial temperature profile
• The solid curve shows how temperature evolves over time
• Heat flows from hot to cold regions (diffusion)
• Sharp features smooth out as higher modes decay faster
• Boundaries are held at u = 0 (like ice baths)

Try These Experiments

• Step function: Watch the Gibbs phenomenon at t=0
• Low modes: See why more terms give better accuracy
• High α: Faster diffusion, quicker smoothing
• Sine wave: Only one mode—pure exponential decay

🤖 Machine Learning Connection

The heat equation is the foundation of diffusion models used in image generation (like Stable Diffusion and DALL-E). These models learn to reverse a diffusion process—starting from noise and "un-diffusing" to create coherent images. The Fourier analysis here explains why diffusion naturally removes high-frequency noise while preserving low-frequency structure.

Machine Learning Connections

The mathematics of Fourier series and the heat equation has profound connections to modern machine learning, particularly in generative AI.

1. Diffusion Models (DALL-E, Stable Diffusion, Midjourney)

Diffusion models work by:

Forward process: Gradually add Gaussian noise to images, simulating heat diffusion that destroys information
Reverse process: Train a neural network to "undo" the diffusion, recovering structure from noise

The Fourier analysis shows why this works: diffusion removes high-frequency details first, leaving low-frequency structure. The network learns to restore frequencies in the reverse order.

2. Score-Based Generative Models

The mathematical foundation of diffusion models is the score function ∇_x log p(x, t), which satisfies a PDE related to the heat equation. The Fokker-Planck equation governing diffusion is:

\frac{\partial p}{\partial t} = \nabla \cdot (p \nabla V) + \sigma^2 \nabla^2 p

The Laplacian term ∇²p is exactly the heat equation! This explains why diffusion gradually "forgets" the original distribution.

3. Spectral Methods in Neural Networks

Spectral normalization: Controls the singular values (eigenvalues) of weight matrices for stable training
Fourier neural operators: Learn in the frequency domain for efficient PDE solving
Graph neural networks: Often based on the graph Laplacian, a discrete analog of ∇²

4. Understanding Neural Network Dynamics

The training dynamics of neural networks can be analyzed using Fourier methods. Key insights include:

Networks learn low-frequency patterns first (similar to how diffusion preserves low frequencies)
High-frequency features require more training time or specific architectural choices
Regularization can be understood as damping high-frequency components

Python Implementation

Fourier Series Solution in Python

🐍fourier_heat_solution.py

Explanation(5)

Code(100)

3Fourier Series Solution Function

This function computes the temperature at any point x and time t by summing Fourier modes. Each mode n contributes sin(nπx/L) times an exponential decay factor.

15Computing Fourier Coefficients

The coefficient b_n depends on the initial condition. For a triangle peaked at L/2, the formula involves sin(nπ/2), which alternates between 1, 0, -1, 0 for n = 1, 2, 3, 4...

33Eigenvalue and Decay

λ_n = (nπ/L)² is the eigenvalue for mode n. The decay factor exp(-αλ_n t) means higher modes decay faster because λ_n grows as n².

41Mode Decay Visualization

This function plots how the amplitude of each mode decreases over time. Mode 1 decays slowly; mode 5 decays 25 times faster!

54Time Constant τ_n

Each mode has a characteristic time τ_n = 1/(αλ_n). After time τ_n, the mode's amplitude has decreased to 1/e ≈ 37% of its initial value.

95 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3from matplotlib.animation import FuncAnimation
4
5def fourier_heat_solution(x, t, L=1.0, alpha=0.01, n_modes=50,
6                          initial_condition='triangle'):
7    """
8    Compute the Fourier series solution to the 1D heat equation.
9
10    ∂u/∂t = α ∂²u/∂x²
11
12    with u(0,t) = u(L,t) = 0 (Dirichlet boundary conditions)
13    """
14    u = np.zeros_like(x)
15
16    for n in range(1, n_modes + 1):
17        # Compute Fourier coefficient b_n
18        if initial_condition == 'triangle':
19            # Triangle peaked at L/2: f(x) = 2x if x < L/2, 2(L-x) otherwise
20            if n % 2 == 0:
21                bn = 0
22            else:
23                bn = (8 / (n**2 * np.pi**2)) * np.sin(n * np.pi / 2)
24        elif initial_condition == 'sine':
25            # sin(πx/L)
26            bn = 1.0 if n == 1 else 0.0
27        elif initial_condition == 'step':
28            # Step function: 1 for 0.25L < x < 0.75L
29            if n % 2 == 0:
30                bn = 0
31            else:
32                npi = n * np.pi
33                bn = (4/npi) * (np.sin(3*npi/4) - np.sin(npi/4))
34        else:
35            bn = 0
36
37        if abs(bn) < 1e-12:
38            continue
39
40        # Eigenvalue and decay factor
41        lambda_n = (n * np.pi / L)**2
42        decay = np.exp(-alpha * lambda_n * t)
43
44        # Add this mode's contribution
45        u += bn * np.sin(n * np.pi * x / L) * decay
46
47    return u
48
49
50def visualize_modal_decay(L=1.0, alpha=0.01, n_modes=5, t_max=2.0):
51    """
52    Visualize how individual Fourier modes decay over time.
53    """
54    t_values = np.linspace(0, t_max, 100)
55
56    plt.figure(figsize=(10, 5))
57
58    for n in range(1, n_modes + 1):
59        lambda_n = (n * np.pi / L)**2
60        decay = np.exp(-alpha * lambda_n * t_values)
61        plt.plot(t_values, decay, label=f'Mode n={n}', linewidth=2)
62
63    plt.xlabel('Time t')
64    plt.ylabel('Amplitude decay factor')
65    plt.title('Exponential Decay of Fourier Modes')
66    plt.legend()
67    plt.grid(True, alpha=0.3)
68    plt.show()
69
70    # Print decay time constants
71    print("\nMode decay time constants τ_n = 1/(αλ_n):")
72    for n in range(1, n_modes + 1):
73        lambda_n = (n * np.pi / L)**2
74        tau = 1 / (alpha * lambda_n)
75        print(f"  n={n}: τ_{n} = {tau:.4f}, λ_{n} = {lambda_n:.4f}")
76
77
78# Demonstration
79L = 1.0
80x = np.linspace(0, L, 200)
81
82# Compare different times
83times = [0, 0.1, 0.5, 1.0, 2.0]
84plt.figure(figsize=(12, 5))
85
86for t in times:
87    u = fourier_heat_solution(x, t, L=L, alpha=0.02, n_modes=50,
88                               initial_condition='triangle')
89    alpha_val = 0.3 + 0.7 * (1 - t/max(times)) if t > 0 else 1.0
90    plt.plot(x, u, label=f't = {t}', alpha=alpha_val, linewidth=2)
91
92plt.xlabel('Position x')
93plt.ylabel('Temperature u(x,t)')
94plt.title('Heat Equation Solution: Fourier Series with 50 Modes')
95plt.legend()
96plt.grid(True, alpha=0.3)
97plt.show()
98
99# Visualize mode decay
100visualize_modal_decay(n_modes=5, t_max=2.0)

Common Pitfalls

Truncating the Series Too Early

Using too few Fourier terms can lead to poor accuracy, especially for initial conditions with sharp features (like step functions). The Gibbs phenomenon causes ~9% overshoot near discontinuities regardless of how many terms you use.

Forgetting the Boundary Conditions

The sine series is appropriate for Dirichlet (fixed value) boundary conditions at both ends. For Neumann (fixed flux) or mixed conditions, you need different eigenfunctions (cosines or combinations).

Confusing Eigenvalues and Coefficients

The eigenvalues λ_n = (nπ/L)² determine how fasteach mode decays. The coefficients b_n determine how much of each mode is initially present. Don't confuse them!

Numerical Verification

Always verify your Fourier series solution against:

The initial condition at t = 0 (sum should equal f(x))
The boundary conditions at all times
Conservation of energy (integral of u² decreases monotonically)
Finite difference solutions for comparison

Test Your Understanding

📝Test Your UnderstandingQuestion 1 of 8

Why do higher Fourier modes decay faster in the heat equation solution?

Summary

Fourier series solutions reveal the deep structure of the heat equation. By decomposing the initial temperature into eigenmodes, we transform a complex PDE into an infinite collection of simple exponential decays.

Key Concepts

Concept	Description
Separation of variables	Assume u(x,t) = X(x)·T(t) to decouple the PDE
Eigenvalue problem	X'' = -λX with BCs → discrete eigenvalues λₙ = (nπ/L)²
Eigenfunctions	Xₙ(x) = sin(nπx/L) satisfy BCs and form a complete basis
Modal decay	Each mode decays as exp(-αλₙt); higher modes decay faster
Fourier coefficients	bₙ = (2/L)∫f(x)sin(nπx/L)dx projects initial condition
Superposition	General solution is sum of all modes: Σ bₙXₙ(x)Tₙ(t)

Key Takeaways

Separation of variables reduces the heat equation PDE to two ODEs: one in space (eigenvalue problem) and one in time (exponential decay)
The boundary conditions quantize the allowed solutions, giving discrete eigenvalues λₙ = (nπ/L)²
Higher modes decay quadratically faster: mode n decays n² times faster than mode 1
This differential decay explains why the heat equation smooths temperature profiles—sharp features disappear first
The same mathematics underlies diffusion models in generative AI, explaining why noise is added/removed in a specific order
The Fourier coefficients completely determine the solution for all future times—all information is in the initial condition

The Essence of Fourier Solutions:

"Any initial temperature distribution can be written as a sum of sine waves. Each wave decays exponentially at its own rate. High frequencies fade fast; low frequencies persist."

Coming Next: In the next section, we'll explore Steady-State Solutions—what happens as t → ∞? When the transient modes have all decayed, what remains is the equilibrium temperature distribution, governed by Laplace's equation.

Learning Objectives

The Big Picture: Why Fourier Series?

The Core Strategy

Historical Context

1807: Fourier's Revolution

The Physical Insight

Modern Legacy

Separation of Variables

Step 1: Start with the Heat Equation

Key Equation

💡 Key Insight

Summary of the Method

Spatial Part: Xn(x)

Temporal Part: Tn(t)

The Eigenvalue Problem

The Eigenvalue Problem

Why Discrete Eigenvalues?

Orthogonality

The Fourier Series Solution

The Fourier Series Solution

What This Solution Tells Us

Computing Fourier Coefficients

Fourier Coefficient Formula

Original Function vs Fourier Approximation

Fourier Sine Coefficients bn

The Fourier Sine Coefficients

Computing Coefficients

Key Insights

Example: Triangular Initial Condition

Pattern in Coefficients

Modal Decay and Physical Interpretation

Mode Decay Rates

Key Observations

The Fourier Series Solution

Why Do Higher Modes Decay Faster?

Mathematical Reason

Physical Reason

The Smoothing Effect

Interactive Solution Explorer

Space-Time Evolution (Heatmap)

What You're Seeing

Try These Experiments

🤖 Machine Learning Connection

Machine Learning Connections

1. Diffusion Models (DALL-E, Stable Diffusion, Midjourney)

2. Score-Based Generative Models

3. Spectral Methods in Neural Networks

4. Understanding Neural Network Dynamics

Python Implementation

Common Pitfalls

Truncating the Series Too Early

Forgetting the Boundary Conditions

Confusing Eigenvalues and Coefficients

Numerical Verification

Test Your Understanding

Summary

Key Concepts

Key Takeaways

Spatial Part: X_n(x)

Temporal Part: T_n(t)

Fourier Sine Coefficients b_n