Solver parameters
Formula and discretization
The heat equation (diffusion) in 1D:
u_t = α u_{xx}, for x ∈ (0,L), t > 0
Explicit (FTCS) update:
u_i^{n+1} = u_i^n + r (u_{i+1}^n - 2 u_i^n + u_{i-1}^n), r = α Δt / Δx^2
Crank–Nicolson (implicit, second order):
(I - θ r A) u^{n+1} = (I + (1-θ) r A) u^n, θ = 0.5 for CN
Crank–Nicolson is unconditionally stable and more accurate in time.
How to use
- Set domain length, material diffusivity α, grid resolution Nx and time steps Nt.
- Choose a numerical scheme: FTCS explicit or Crank–Nicolson implicit.
- Select initial condition and boundary conditions; adjust parameters if needed.
- Press Solve. Results appear above the form with plots and sample tables.
- Use Download CSV to export data or Export PDF to print a report.
Example data
Default example uses L=1.0, α=0.01, Nx=50, T=0.5, Nt=500, sinusoidal initial condition u(x,0)=sin(πx/L).
| x | u(x,0) |
|---|---|
| 0 | 0 |
| 0.1 | 0.309017 |
| 0.2 | 0.587785 |
| 0.3 | 0.809017 |
| 0.4 | 0.951057 |
| 0.5 | 1 |
| 0.6 | 0.951057 |
| 0.7 | 0.809017 |
| 0.8 | 0.587785 |
| 0.9 | 0.309017 |
| 1 | 0 |
Understanding transient heat diffusion
The transient heat equation models how temperature evolves due to diffusion. In one dimension it is u_t = α u_{xx}; α (m²/s) controls rate. For practical engineering, typical α values range from 1e-7 (insulators) to 1e-4 (metals). The solver supports explicit and Crank–Nicolson schemes to balance accuracy and stability.
Numerical stability and timestep selection
Explicit FTCS requires r = α Δt / Δx² ≤ 0.5 for stability. For example, with L=1.0, Nx=50, α=0.01, Δx=0.02, the explicit limit gives Δt ≤ 0.5*(Δx²)/α ≈ 0.02. The Crank–Nicolson scheme is unconditionally stable and allows larger time steps while maintaining second‑order temporal accuracy.
Discretization and error
Spatial discretization error is O(Δx²) for central differences. Time discretization error is O(Δt) for explicit Euler and O(Δt²) for Crank–Nicolson. Reduce Δx and Δt to improve accuracy; monitor convergence by halving Δx and comparing norms of the numerical solution.
Boundary and initial conditions
Dirichlet boundaries fix temperature values at domain ends; Neumann boundaries model insulated ends (zero flux). Initial conditions control transient behavior — common choices include sinusoidal profiles (useful for modal tests) and Gaussian pulses (useful for localized heat sources).
Practical example
Using the built‑in example: L=1.0, α=0.01, Nx=50, Nt=500, T=0.5 with u(x,0)=sin(πx) produces smooth decay of the fundamental mode. Snapshots at evenly spaced times show amplitude reduction consistent with analytical modal decay rates.
Performance considerations
Interior tridiagonal solves in Crank–Nicolson cost O(Nx) per time step and are memory efficient. For large grids, consider compiled solvers or multi‑threading. For 2D extensions, matrix size and solver choice substantially affect runtime and memory usage.
Exporting and post‑processing
Use the CSV export to post‑process in Python, MATLAB, or Excel. The PDF export captures plots for reporting. To compute energy or L2 norms, integrate u(x,t)² over the domain using simple trapezoidal sums from exported data.
When to use each scheme
Choose FTCS for quick tests with fine temporal resolution and small Δt. Choose Crank–Nicolson for production simulations requiring larger Δt and improved accuracy. Combine grid refinement and time‑step studies to validate results.
Frequently asked questions
1. Is Crank–Nicolson always better than explicit FTCS?
Crank–Nicolson is unconditionally stable and more accurate in time, but requires solving a linear system each step. FTCS is simpler and cheaper per step but requires small Δt for stability.
2. How do I choose Nx and Nt?
Balance resolution and cost: halve Δx to test spatial convergence; choose Δt to satisfy stability for FTCS or to match desired temporal resolution for implicit schemes.
3. Can this solver handle variable diffusivity?
The current version assumes constant α. Variable diffusivity requires modifying coefficients in the discretization and updating the tridiagonal assembly accordingly.
4. How accurate is the Gaussian initial condition?
A Gaussian pulse is smooth and well‑resolved with moderate grids. Accuracy depends on σ relative to Δx; ensure at least 6–10 points across the Gaussian width for good resolution.
5. Can I extend to 2D or 3D?
Yes. Extension requires 2D finite difference stencils and solvers (e.g., ADI or iterative Krylov methods). Memory and compute requirements grow rapidly with dimension.
6. What if I need Neumann boundaries?
The solver includes a simple zero‑flux Neumann approximation. For higher accuracy, implement ghost points or second‑order Neumann discretizations.
7. How do I verify correctness?
Compare against analytical solutions (e.g., sinusoidal modes) or perform grid/time refinement studies and check convergence rates consistent with theoretical order.