Numerical simulation of heat diffusion in a 1D bar with the FTCS scheme
Papadopol Lucian Ioan
l.i.papadopol@gmail.com
1. Objective of the simulation
In this exercise we simulate heat diffusion in a one-dimensional metal bar.
We imagine a thin bar, initially cold, in which a central region is heated. After the initial heating, heat tends to propagate from the hotter regions toward the colder regions.
The goal of the simulation is to observe how the temperature profile along the bar changes over time.
The main physical phenomenon is diffusion, that is, the process through which a quantity progressively spreads in space.
In this case, the quantity that diffuses is temperature.
2. Physical system
We consider a thin metal bar of length L.
The bar is described as a one-dimensional domain:
0 \le x \le L
The unknown quantity is:
u(x,t)
where:
u(x,t) represents the temperature at position x and time t;
x is the position along the bar;
t is time.
We assume that initially the bar is cold, except for a hotter central region.
As time passes, heat diffuses toward the colder regions.
Physically, we expect that:
the temperature peak decreases;
the hot region widens;
the profile becomes progressively smoother;
temperature distributes along the bar.
3. Difference between convection and diffusion
In the previous notebooks we studied convection.
For linear convection, the model was:
\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0
In that case, the profile was transported along the domain.
For nonlinear convection, the model was:
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=0
In that case, the profile was transported and deformed, because the propagation speed depended on the local value of u.
Now the phenomenon is different.
In diffusion, the profile is not simply transported to the right or to the left.
Instead, the profile tends to spread out in space.
Diffusion is described by an equation containing a second derivative with respect to space.
4. The 1D diffusion equation
The mathematical model for heat diffusion in a 1D bar is:
\frac{\partial u}{\partial t}=\alpha\frac{\partial^2 u}{\partial x^2}
This is the one-dimensional heat equation.
The parameter:
\alpha
is the thermal diffusivity.
Thermal diffusivity measures how rapidly heat diffuses in the material.
If \alpha is large, heat spreads rapidly. If \alpha is small, heat diffuses more slowly.
5. Interpretation of the terms
The equation is:
\frac{\partial u}{\partial t}=\alpha\frac{\partial^2 u}{\partial x^2}
The term:
\frac{\partial u}{\partial t}
describes the time variation of temperature.
The term:
\frac{\partial^2 u}{\partial x^2}
describes the spatial curvature of the temperature profile.
Diffusion is directly related to curvature.
If the profile has strong curvature, this means that there are large temperature differences between neighboring points. In that case, heat tends to redistribute rapidly.
If the profile is almost flat, the curvature is small and the temperature changes slowly.
6. Spatial domain
We consider a bar of length L.
The spatial domain is:
0 \le x \le L
We divide the bar into a grid of points:
x_i=i\Delta x
where:
i=0,1,2,\dots,N-1
The spatial step is:
\Delta x=\frac{L}{N-1}
The numerical value:
u_i^n
approximates:
u(x_i,t_n)
that is, the temperature at grid point x_i and time t_n.
7. Time discretization
Time is also discretized:
t_n=n\Delta t
where:
n=0,1,2,\dots
and \Delta t is the time step.
The numerical solution is therefore computed on a space-time grid.
Each time row contains the temperature profile along the bar at one instant:
time t0 -> profile u(x,t0)
time t1 -> profile u(x,t1)
time t2 -> profile u(x,t2)
...
8. FTCS scheme
To solve the diffusion equation numerically, we use the FTCS scheme.
FTCS means:
Forward Time;
Central Space.
We therefore use:
a forward difference for the time derivative;
a centered difference for the second spatial derivative.
The starting equation is:
\frac{\partial u}{\partial t}=\alpha\frac{\partial^2 u}{\partial x^2}
The time derivative is approximated as:
\frac{\partial u}{\partial t}\approx\frac{u_i^{n+1}-u_i^n}{\Delta t}
The second spatial derivative is approximated with a centered difference:
\frac{\partial^2 u}{\partial x^2}\approx\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\Delta x^2}
Substituting into the diffusion equation gives:
\frac{u_i^{n+1}-u_i^n}{\Delta t}=\alpha\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\Delta x^2}
Solving for u_i^{n+1} gives:
u_i^{n+1}=u_i^n+\alpha\frac{\Delta t}{\Delta x^2}\left(u_{i+1}^n-2u_i^n+u_{i-1}^n\right)
This is the explicit FTCS scheme for 1D diffusion.
9. Diffusive number
We introduce the parameter:
r=\alpha\frac{\Delta t}{\Delta x^2}
This parameter measures the relation between the time step, the spatial step and the material diffusivity.
With this definition, the FTCS scheme becomes:
u_i^{n+1}=u_i^n+r\left(u_{i+1}^n-2u_i^n+u_{i-1}^n\right)
or equivalently:
u_i^{n+1}=(1-2r)u_i^n+r u_{i+1}^n+r u_{i-1}^n
This form clearly shows that the new value at point i depends on the current value and on the two neighboring points.
10. Stability condition
The FTCS scheme for the 1D diffusion equation is explicit.
This means that the profile at the next time level is computed directly from the current profile.
However, as in the convection case, the time step cannot be chosen arbitrarily.
For 1D diffusion with the FTCS scheme, the stability condition is:
r\le\frac{1}{2}
that is:
\alpha\frac{\Delta t}{\Delta x^2}\le\frac{1}{2}
From this relation we obtain the time-step limit:
\Delta t\le\frac{\Delta x^2}{2\alpha}
This condition is very important. If the time step is too large, the simulation can become unstable and produce nonphysical oscillations.
11. Difference between convective CFL and diffusive stability
In the convection case, stability was controlled by the Courant number:
C=\frac{|c|\Delta t}{\Delta x}
or, in the nonlinear case:
C=\frac{\max |u|\Delta t}{\Delta x}
In the diffusion case, the stability parameter is instead:
r=\alpha\frac{\Delta t}{\Delta x^2}
The fundamental difference is that convection involves \Delta x, while diffusion involves \Delta x^2.
This means that, if the spatial grid is refined, the time step for diffusion must decrease very rapidly.
12. Boundary conditions
In this simulation we consider a bar with endpoints kept cold.
We impose:
u(0,t)=0
and:
u(L,t)=0
These are Dirichlet boundary conditions.
Physically, they mean that the endpoints of the bar are kept at zero temperature, for example because they are connected to two cold thermal reservoirs.
In the code we impose:
u_0^n=0
and:
u_{N-1}^n=0
for every time level.
13. Initial condition
The initial condition describes the temperature of the bar at:
t=0
We use a localized profile in the central part of the bar.
We can imagine that a central region is heated, while the rest of the bar remains cold.
A convenient choice is a Gaussian bell:
u(x,0)=\exp\left(-\frac{(x-x_c)^2}{2\sigma^2}\right)
where:
x_c is the center of the hot region;
\sigma controls the width of the heated region.
This choice is useful because it produces a smooth profile.
During the simulation, we expect the peak to decrease and the bell to widen.
14. What we expect from the simulation
From the simulation we expect that:
the initial hot region widens over time;
the maximum temperature decreases;
the profile becomes smoother;
heat diffuses toward the colder regions;
the endpoints remain cold because of the boundary conditions.
Unlike convection, the profile should not simply move.
In the diffusive case, the main phenomenon is the smoothing of the profile.
15. Numerical objective of the notebook
In the notebook we want to:
build the spatial grid;
define a localized initial temperature;
choose a time step compatible with the stability condition;
implement the explicit FTCS scheme;
impose the boundary conditions;
observe the evolution of the temperature profile;
verify that the profile widens and decreases over time.
The main result will be the graph of:
u(x,t)
at different times.
This allows us to observe the typical behavior of diffusive, or parabolic, PDEs.
First step — Initial conditions
import numpy as np
import matplotlib.pyplot as plt
# Bar length
L = 1.0
# Number of grid points
nx = 201
# Spatial grid
x = np.linspace(0.0, L, nx)
# Spatial step
dx = x[1] - x[0]
# Thermal diffusivity
alpha = 0.01
# Desired diffusive number
r_target = 0.4
# Time step chosen to satisfy the FTCS stability condition
dt = r_target * dx**2 / alpha
# Effective diffusive number
r = alpha * dt / dx**2
print("dx =", dx)
print("dt =", dt)
print("r =", r)
# Initial condition: cold bar with a hot central region
x_c = 0.5
sigma = 0.06
u0 = np.exp(-((x - x_c)**2) / (2 * sigma**2))
# Boundary conditions: cold endpoints
u0[0] = 0.0
u0[-1] = 0.0
plt.figure(figsize=(8, 4))
plt.plot(x, u0, label="Initial temperature")
plt.xlabel("x")
plt.ylabel("u(x,0)")
plt.title("Initial temperature in the bar")
plt.grid(True)
plt.legend()
plt.show()
Initial temperature in the bar.
The Gaussian bell represents a central hot region in an otherwise cold bar. The boundary values are set to zero, meaning that the endpoints are kept cold.
Second step — Definition of the simulation function
def diffusion_ftcs(u0, alpha, dx, dt, nt, left_bc=None, right_bc=None):
u = np.asarray(u0, dtype=float).copy()
r = alpha * dt / dx**2
if r > 0.5:
raise ValueError("Unstable FTCS choice")
for _ in range(nt):
un = u.copy()
u[1:-1] = un[1:-1] + r * (
un[2:] - 2 * un[1:-1] + un[:-2]
)
if left_bc is not None:
u[0] = left_bc
if right_bc is not None:
u[-1] = right_bc
return u
The central update line implements the FTCS formula:
u_i^{n+1}=u_i^n+r\left(u_{i+1}^n-2u_i^n+u_{i-1}^n\right)
The stability check prevents running the explicit method with r>0.5.
Third step — Simulation
# Final simulation time
t_end = 1.0
# Number of time steps
nt = int(t_end / dt)
print("nt =", nt)
print("simulated final time =", nt * dt)
# Simulation
u_final = diffusion_ftcs(
u0,
alpha,
dx,
dt,
nt,
left_bc=0.0,
right_bc=0.0
)
With the parameters used here, the numerical values are:
dx = 0.005
dt = 0.001
r = 0.4
nt = 1000
simulated final time = 1.0
Initial and final temperature profiles.
The central bell:
decreases in height;
widens;
becomes smoother;
remains zero at the endpoints.
This is exactly the expected behavior of diffusion. The heat does not travel as a rigid profile; instead, it spreads out from the hotter region toward the colder parts of the bar.
