36_heat_equation_1d/
36_heat_equation_1d.rs

1//! # Example: 1D Heat Equation
2//!
3//! Run: cargo run --example 36_heat_equation_1d
4//!
5//! ## Problem
6//! A thin rod has its two ends held at fixed temperatures and starts cold inside.
7//! How does the temperature profile evolve, and where does it settle?
8//!
9//! ## Math idea
10//! The explicit (forward-Euler) finite-difference update of the 1D heat equation is
11//! `u_i ← u_i + r·(u_{i-1} - 2·u_i + u_{i+1})`. Collecting one full sweep into a
12//! tridiagonal matrix `A` (with identity rows at the fixed boundaries) makes each
13//! time step a single `u_next = A · u`. The stability condition is `r ≤ 0.5`.
14//!
15//! ## Tensor representation
16//! `A` is an N×N `Tensor`; the rod state `u` is length-N. Each time step is one
17//! `Tensor::matmul` (`[n, n] × [n] -> [n]`).
18//!
19//! ## What this demonstrates
20//! - encoding a stencil as a matrix and iterating with `Tensor::matmul`;
21//! - fixed boundary conditions via identity rows;
22//! - convergence to the steady-state (linear) profile between the boundaries.
23//!
24//! ## Expected output
25//! ```text
26//! step    u(rod)
27//!    1   [0.00, 0.00, 0.00, 25.00, 100.00]
28//!   10   [0.00, 16.26, 37.61, 66.22, 100.00]
29//!   50   [0.00, 24.98, 49.98, 74.98, 100.00]
30//!  300   [0.00, 25.00, 50.00, 75.00, 100.00]
31//! Heat equation 1D: OK
32//! ```
33
34use matten::Tensor;
35
36/// One explicit heat-equation step: `u_next = A · u`.
37fn step(a: &Tensor, u: &Tensor) -> Tensor {
38    a.matmul(u)
39}
40
41fn main() {
42    // 5 nodes along the rod. Boundaries are fixed (left = 0, right = 100); the
43    // interior starts cold. Diffusion number r = 0.25 satisfies r <= 0.5 (stable).
44    let r = 0.25;
45    let d = 1.0 - 2.0 * r;
46
47    // Boundary rows are identity (hold the value); interior rows are [r, 1-2r, r].
48    let a = Tensor::new(
49        vec![
50            1.0, 0.0, 0.0, 0.0, 0.0, //
51            r, d, r, 0.0, 0.0, //
52            0.0, r, d, r, 0.0, //
53            0.0, 0.0, r, d, r, //
54            0.0, 0.0, 0.0, 0.0, 1.0, //
55        ],
56        &[5, 5],
57    );
58
59    let mut u = Tensor::from_vec(vec![0.0, 0.0, 0.0, 0.0, 100.0]);
60
61    println!("step    u(rod)");
62    for s in 1..=300 {
63        u = step(&a, &u);
64        if matches!(s, 1 | 10 | 50 | 300) {
65            let row: Vec<String> = u.as_slice().iter().map(|v| format!("{v:.2}")).collect();
66            println!("{s:>4}   [{}]", row.join(", "));
67        }
68    }
69
70    // Steady state is the linear profile between the two boundary temperatures.
71    let expected = [0.0, 25.0, 50.0, 75.0, 100.0];
72    for (got, want) in u.as_slice().iter().zip(&expected) {
73        assert!((got - want).abs() < 0.1);
74    }
75    println!("Heat equation 1D: OK");
76}
36_heat_equation_1d/36_heat_equation_1d.rs

36_heat_equation_1d/
36_heat_equation_1d.rs
