Tutorial: Solution of the heat equation with Dirichlet boundary conditions
We continue the previous tutorial on solving the heat equation with Neumann boundary conditions by looking at Dirichlet boundary conditions instead, resulting in a non-conservative production-destruction system.
Definition of the (non-conservative) production-destruction system
Consider the heat equation
\[\partial_t u(t,x) = \mu \partial_x^2 u(t,x),\quad u(0,x)=u_0(x),\]
with $μ ≥ 0$, $t≥ 0$, $x\in[0,1]$, and homogeneous Dirichlet boundary conditions. We use again a finite volume discretization, i.e., we split the domain $[0, 1]$ into $N$ uniform cells of width $\Delta x = 1 / N$. As degrees of freedom, we use the mean values of $u(t)$ in each cell approximated by the point value $u_i(t)$ in the center of cell $i$. Finally, we use the classical central finite difference discretization of the Laplacian with homogeneous Dirichlet boundary conditions, resulting in the ODE
\[\partial_t u(t) = L u(t), \quad L = \frac{\mu}{\Delta x^2} \begin{pmatrix} -2 & 1 \\ 1 & -2 & 1 \\ & \ddots & \ddots & \ddots \\ && 1 & -2 & 1 \\ &&& 1 & -2 \end{pmatrix}.\]
The system can be written as a non-conservative PDS with production terms
\[\begin{aligned} &p_{i,i-1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i-1}(t),\quad i=2,\dots,N, \\ &p_{i,i+1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i+1}(t),\quad i=1,\dots,N-1, \end{aligned}\]
and destruction terms $d_{i,j} = p_{j,i}$ for $i \ne j$ as well as the non-conservative destruction terms
\[\begin{aligned} d_{1,1}(t,\mathbf u(t)) &= \frac{\mu}{\Delta x^2} u_{1}(t), \\ d_{N,N}(t,\mathbf u(t)) &= \frac{\mu}{\Delta x^2} u_{N}(t). \end{aligned}\]
In addition, all production and destruction terms not listed are zero.
Solution of the non-conservative production-destruction system
Now we are ready to define a PDSProblem and to solve this problem with a method of PositiveIntegrators.jl or OrdinaryDiffEq.jl. In the following we use $N = 100$ nodes and the time domain $t \in [0,1]$. Moreover, we choose the initial condition
\[u_0(x) = \sin(\pi x)^2.\]
x_boundaries = range(0, 1, length = 101)
x = x_boundaries[1:end-1] .+ step(x_boundaries) / 2
u0 = @. sinpi(x)^2 # initial solution
tspan = (0.0, 1.0) # time domainWe will choose three different matrix types for the production terms and the resulting linear systems:
- standard dense matrices (default)
- sparse matrices (from SparseArrays.jl)
- tridiagonal matrices (from LinearAlgebra.jl)
Standard dense matrices
using PositiveIntegrators # load ConservativePDSProblem
function heat_eq_P!(P, u, μ, t)
fill!(P, 0)
N = length(u)
Δx = 1 / N
μ_Δx2 = μ / Δx^2
let i = 1
# Dirichlet boundary condition
P[i, i + 1] = u[i + 1] * μ_Δx2
end
for i in 2:(length(u) - 1)
# interior stencil
P[i, i - 1] = u[i - 1] * μ_Δx2
P[i, i + 1] = u[i + 1] * μ_Δx2
end
let i = length(u)
# Dirichlet boundary condition
P[i, i - 1] = u[i - 1] * μ_Δx2
end
return nothing
end
function heat_eq_D!(D, u, μ, t)
fill!(D, 0)
N = length(u)
Δx = 1 / N
μ_Δx2 = μ / Δx^2
# Dirichlet boundary condition
D[begin] = u[begin] * μ_Δx2
D[end] = u[end] * μ_Δx2
return nothing
end
μ = 1.0e-2
prob = PDSProblem(heat_eq_P!, heat_eq_D!, u0, tspan, μ) # create the PDS
sol = solve(prob, MPRK22(1.0); save_everystep = false)using Plots
plot(x, u0; label = "u0", xguide = "x", yguide = "u")
plot!(x, last(sol.u); label = "u")
Sparse matrices
To use different matrix types for the production terms and linear systems, you can use the keyword argument p_prototype of ConservativePDSProblem and PDSProblem.
using SparseArrays
p_prototype = spdiagm(-1 => ones(eltype(u0), length(u0) - 1),
+1 => ones(eltype(u0), length(u0) - 1))
prob_sparse = PDSProblem(heat_eq_P!, heat_eq_D!, u0, tspan, μ;
p_prototype = p_prototype)
sol_sparse = solve(prob_sparse, MPRK22(1.0); save_everystep = false)plot(x, u0; label = "u0", xguide = "x", yguide = "u")
plot!(x, last(sol_sparse.u); label = "u")
Tridiagonal matrices
The sparse matrices used in this case have a very special structure since they are in fact tridiagonal matrices. Thus, we can also use the special matrix type Tridiagonal from the standard library LinearAlgebra.
using LinearAlgebra
p_prototype = Tridiagonal(ones(eltype(u0), length(u0) - 1),
ones(eltype(u0), length(u0)),
ones(eltype(u0), length(u0) - 1))
prob_tridiagonal = PDSProblem(heat_eq_P!, heat_eq_D!, u0, tspan, μ;
p_prototype = p_prototype)
sol_tridiagonal = solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false)plot(x, u0; label = "u0", xguide = "x", yguide = "u")
plot!(x, last(sol_tridiagonal.u); label = "u")
Performance comparison
Finally, we use BenchmarkTools.jl to compare the performance of the different implementations.
using BenchmarkTools
@benchmark solve(prob, MPRK22(1.0); save_everystep = false)BenchmarkTools.Trial: 534 samples with 1 evaluation per sample.
Range (min … max): 8.669 ms … 16.817 ms ┊ GC (min … max): 0.00% … 2.86%
Time (median): 9.086 ms ┊ GC (median): 3.75%
Time (mean ± σ): 9.353 ms ± 986.971 μs ┊ GC (mean ± σ): 3.61% ± 3.02%
▂ ▂█▆▄▂▃▁
▇██████████▆▆▅▅▆▁▄▁▅▁▄▁▄▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▄▄▁▁▄▁▄▄▁▄▁▄▁▄▇▆▁▅▁▅ ▇
8.67 ms Histogram: log(frequency) by time 13.7 ms <
Memory estimate: 10.21 MiB, allocs estimate: 712.@benchmark solve(prob_sparse, MPRK22(1.0); save_everystep = false)BenchmarkTools.Trial: 746 samples with 1 evaluation per sample.
Range (min … max): 6.215 ms … 16.890 ms ┊ GC (min … max): 0.00% … 2.63%
Time (median): 6.762 ms ┊ GC (median): 5.04%
Time (mean ± σ): 6.704 ms ± 714.901 μs ┊ GC (mean ± σ): 3.50% ± 2.57%
▆▄▂▃ ▅█▃
████▇▅▄▅███▆▇▄▅▄▄▅▄▁▄▁▁▁▁▁▁▅▁▄▁▄▁▁▁▄▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▁▁▄ ▇
6.22 ms Histogram: log(frequency) by time 9.78 ms <
Memory estimate: 10.19 MiB, allocs estimate: 5338.By default, we use an LU factorization for the linear systems. At the time of writing, Julia uses SparseArrays.jl defaulting to UMFPACK from SuiteSparse in this case. However, the linear systems do not necessarily have the structure for which UMFPACK is optimized for. Thus, it is often possible to gain performance by switching to KLU instead.
using LinearSolve
@benchmark solve(prob_sparse, MPRK22(1.0; linsolve = KLUFactorization()); save_everystep = false)BenchmarkTools.Trial: 3540 samples with 1 evaluation per sample.
Range (min … max): 1.267 ms … 14.540 ms ┊ GC (min … max): 0.00% … 10.90%
Time (median): 1.350 ms ┊ GC (median): 0.00%
Time (mean ± σ): 1.410 ms ± 438.009 μs ┊ GC (mean ± σ): 2.10% ± 5.79%
▁▄▆██▇▅▂▁ ▄▃ ▁
██████████▇▇██▇▅▆▅▆▃▃▅▃▄▃▁▁▁▃▁▁▁▁▄▃▁▁▃▃▄▁▆▃▃▃▁▄▁▅▅▆▇▆▇▆▇▆▃▅ █
1.27 ms Histogram: log(frequency) by time 2.35 ms <
Memory estimate: 542.10 KiB, allocs estimate: 1078.@benchmark solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false)BenchmarkTools.Trial: 9336 samples with 1 evaluation per sample.
Range (min … max): 472.752 μs … 23.595 ms ┊ GC (min … max): 0.00% … 97.80%
Time (median): 487.059 μs ┊ GC (median): 0.00%
Time (mean ± σ): 533.275 μs ± 408.158 μs ┊ GC (mean ± σ): 4.48% ± 8.78%
██▇▆▅▃▁ ▂▄▄▂ ▂▂▁ ▂
███████▇▆▇▅▅▄▃▃▄▁▁▄▃▁▆▇▅▄▁▄▁▁▁▁█████▇▆▅▁▅▃▁▅▄▁▁▄▃▁▁▁▁▁▁▆█████ █
473 μs Histogram: log(frequency) by time 933 μs <
Memory estimate: 591.77 KiB, allocs estimate: 1630.Package versions
These results were obtained using the following versions.
using InteractiveUtils
versioninfo()
println()
using Pkg
Pkg.status(["PositiveIntegrators", "SparseArrays", "KLU", "LinearSolve", "OrdinaryDiffEq"],
mode=PKGMODE_MANIFEST)Julia Version 1.11.3
Commit d63adeda50d (2025-01-21 19:42 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × AMD EPYC 7763 64-Core Processor
WORD_SIZE: 64
LLVM: libLLVM-16.0.6 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
Status `~/work/PositiveIntegrators.jl/PositiveIntegrators.jl/docs/Manifest.toml`
[ef3ab10e] KLU v0.6.0
[7ed4a6bd] LinearSolve v2.38.0
[1dea7af3] OrdinaryDiffEq v6.90.1
[d1b20bf0] PositiveIntegrators v0.2.6 `~/work/PositiveIntegrators.jl/PositiveIntegrators.jl`
[2f01184e] SparseArrays v1.11.0