┌─────────────────────────────────────────────────────────────┐
│ INVERSE PROBLEM WORKFLOW │
└─────────────────────────────────────────────────────────────┘
1. SETUP
└─> Define domain [0, L], time [0, T]
└─> Create spatial grid (nx points) and time grid (nt steps)
└─> Define initial condition u(x, 0)
└─> Define boundary conditions u(0,t), u(L,t)
2. DIRECT PROBLEM (Generate "True" Data)
└─> Set k_true = 0.01
└─> Solve heat equation: ∂u/∂t = k·∂²u/∂x²
└─> Get full solution u_true(x, t)
└─> Extract measurements at sensor positions
3. INVERSE PROBLEM (Pretend we don't know k)
└─> Initialize k_guess = 0.02
└─> For each candidate k:
├─> Solve forward problem
├─> Compare with measurements
└─> Compute error J(k)
└─> Optimize to find k that minimizes J(k)
4. VERIFICATION
└─> Compare k_estimated with k_true
└─> Solve forward problem with k_estimated
└─> Visualize differences
5. RESULTS
└─> Report k_estimated, error, convergence
└─> Generate 6 diagnostic plots
HeatEquationSolver
│
├── __init__(L, T, nx, nt)
│ └─> Creates spatial grid x and temporal grid t
│
├── solve_direct(k, u_initial, u_left, u_right)
│ ├─> Checks CFL stability condition
│ ├─> Initializes solution array u[nt, nx]
│ ├─> Time stepping loop (explicit FTCS)
│ └─> Returns u(x, t)
│
└── extract_measurements(u, measurement_points)
└─> Returns u[:, measurement_points]
InverseProblemSolver
│
├── __init__(solver, measurement_points, u_measured, ...)
│ └─> Stores problem data and boundary conditions
│
├── objective_function(k)
│ ├─> Solves direct problem with candidate k
│ ├─> Extracts computed values at measurement points
│ ├─> Computes least squares error
│ └─> Returns J(k) = Σ(u_computed - u_measured)²
│
└── solve(k_initial_guess, method)
├─> Calls scipy.optimize.minimize
├─> Tracks optimization history
└─> Returns {'k_estimated', 'objective_value', 'success', ...}
INPUT PROCESSING OUTPUT
───── ────────── ──────
Initial Condition ─┐
│
Boundary Cond. ───┼──> HeatEquationSolver ──> u_true(x,t)
│ │
True k ────────────┘ │
↓
Extract Measurements ──> u_measured
│
↓
k_initial_guess ──┐ │
│ │
u_measured ───────┼──> InverseProblemSolver
│ │
Boundary Cond. ───┘ │
↓
Optimization Loop
┌────────┴────────┐
↓ ↓
Try candidate k Compute J(k)
↓ ↓
Solve forward Compare with
problem measurements
│ │
└────────┬────────┘
↓
Minimize J(k) ──> k_estimated
│
↓
Verification & Plots
Heat Equation: ∂u/∂t = k·∂²u/∂x²
Discretization:
x_i = i·Δx, i = 0, 1, ..., nx-1
t_n = n·Δt, n = 0, 1, ..., nt-1
FTCS Scheme:
u[n+1, i] = u[n, i] + r·(u[n, i+1] - 2·u[n, i] + u[n, i-1])
where r = k·Δt/Δx² (mesh ratio)
Stability: r ≤ 0.5 (CFL condition)
Objective Function:
J(k) = Σᵢ₌₁ⁿ Σₜ₌₁ᵐ [u_computed(xᵢ, tₜ; k) - u_measured(xᵢ, tₜ)]²
Optimization:
k* = argmin J(k)
k > 0
Methods Available:
• Nelder-Mead (simplex, derivative-free)
• Powell (conjugate directions)
• BFGS (quasi-Newton)
• L-BFGS-B (bounded quasi-Newton)
• CG (conjugate gradient)
• TNC (truncated Newton)
• SLSQP (sequential least squares)
The plot_results() function creates a 6-panel figure:
┌─────────────────┬─────────────────┬─────────────────┐
│ Plot 1 (3D) │ Plot 2 (3D) │ Plot 3 (3D) │
│ True Solution │ Reconstructed │ Absolute Error │
│ with k_true │ with k_estimate │ |u_true - u_rec||
├─────────────────┼─────────────────┼─────────────────┤
│ Plot 4 (2D) │ Plot 5 (2D) │ Plot 6 (2D) │
│ Snapshots │ Time Series │ Optimization │
│ at different │ at measurement │ History │
│ times │ points │ J(k) vs k │
└─────────────────┴─────────────────┴─────────────────┘
L = 1.0 # Domain length
T = 0.5 # Final time
nx = 51 # Spatial points
nt = 501 # Time steps
k_true = 0.01 # True thermal diffusivity
n_interior_points = 5 # Interior measurement points
k_initial_guess = 0.02 # Starting guess
method = 'Nelder-Mead' # Optimization method# Grid resolution
nx: 31-101 (higher = more accurate, slower)
nt: 201-1001 (constrained by CFL condition)
# Thermal diffusivity
k: 0.001-0.1 (typical range)
# Measurement points
n_measurements: 5-15 (including boundaries)
# Initial guess
k_guess: 0.5·k_true to 2·k_true (robust)def max_dt(k, dx):
"""Maximum stable time step."""
return 0.5 * dx**2 / k
def min_nt(k, L, nx, T):
"""Minimum time steps for stability."""
dx = L / (nx - 1)
dt_max = max_dt(k, dx)
return int(np.ceil(T / dt_max)) + 1
# Example:
# k=0.01, L=1.0, nx=51, T=0.5
# → min_nt ≈ 251u_initial = lambda x: np.sin(np.pi * x)
u_left = lambda t: 0.0
u_right = lambda t: 0.0
# Result: ~0.0000% erroru_initial = lambda x: np.exp(-((x-0.5)**2) / 0.02)
u_left = lambda t: 0.0
u_right = lambda t: 0.0
# Result: ~0.0001% erroru_initial = lambda x: 1.0 if x < 0.5 else 0.0
u_left = lambda t: 1.0
u_right = lambda t: 0.0
# Result: ~0.01% error (discontinuity challenge)u_initial = lambda x: 0.0
u_left = lambda t: 100 * np.sin(2*np.pi*t)
u_right = lambda t: 0.0
# Result: ~0.0002% error-
Different Initial Conditions (4 tests)
- Sine wave
- Gaussian
- Step function
- Multi-peak
-
Noise Robustness (5 levels)
- 0%, 1%, 3%, 5%, 10% noise
-
Measurement Density (7 configurations)
- 3, 5, 7, 9, 11, 15, 21 points
-
Algorithm Comparison (7 methods)
- Performance metrics
- Convergence analysis
- Timing comparison
-
Initial Guess Sensitivity (7 values)
- From 10× too small to 20× too large
Computational Complexity:
- Direct problem: O(nx · nt) per solve
- Inverse problem: O(n_iter · nx · nt)
- Memory: O(nx · nt)
Typical Performance (Intel i7, 51×501 grid):
- Single forward solve: 5-10 ms
- Inverse problem (20 iterations): 100-200 ms
- Full visualization: 500-800 ms
Scaling:
- 2× grid resolution → 4× time (2D grid)
- 2× optimization iterations → 2× time (linear)
- Alifanov (1994) - Inverse Heat Transfer Problems
- Beck et al. (1985) - Inverse Heat Conduction
- Ozisik & Orlande (2000) - Inverse Heat Transfer
- LeVeque (2007) - Finite Difference Methods for ODEs/PDEs
- Smith (1985) - Numerical Solution of PDEs
- Nocedal & Wright (2006) - Numerical Optimization
- scipy.optimize documentation