2D time-independent heat equation

2D time-independent heat equation#

From the tutorial

Using MethodOfLines.jl (SciML/MethodOfLines.jl) to symbolically define the PDE system and use the finite difference method (FDM) to solve the following PDE:

\[ \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \]

using ModelingToolkit
using MethodOfLines
using DomainSets
using OrdinaryDiffEq
using NonlinearSolve
using Plots

Setup variables and differential operators

@independent_variables x y
@variables u(..)

Dxx = Differential(x)^2
Dyy = Differential(y)^2

Differential(y) ∘ Differential(y)

PDE equation

eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ 0

\[ \begin{equation} \frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} u\left( x, y \right) + \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} u\left( x, y \right) = 0 \end{equation} \]

Boundary conditions

bcs = [u(0, y) ~ x * y,
       u(1, y) ~ x * y,
       u(x, 0) ~ x * y,
       u(x, 1) ~ x * y
]

\[\begin{split} \begin{align} u\left( 0, y \right) &= x y \\ u\left( 1, y \right) &= x y \\ u\left( x, 0 \right) &= x y \\ u\left( x, 1 \right) &= x y \end{align} \end{split}\]

Space and time domains

domains = [ x ∈ Interval(0.0, 1.0), y ∈ Interval(0.0, 1.0)]

2-element Vector{Symbolics.VarDomainPairing}:
 Symbolics.VarDomainPairing(x, 0.0 .. 1.0)
 Symbolics.VarDomainPairing(y, 0.0 .. 1.0)

PDE system#

@named pdesys = PDESystem(eq, bcs, domains, [x, y], [u(x, y)])

\[ \begin{align} \frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} u\left( x, y \right) + \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} u\left( x, y \right) &= 0 \end{align} \]

Discretize the PDE system into an Nonlinear system Pass nothing to the time parameter

@time prob = let dx=0.1
    discretization = MOLFiniteDifference([x=>dx, y=>dx], nothing, approx_order=2, grid_align = MethodOfLines.EdgeAlignedGrid())
    prob = discretize(pdesys, discretization)
end

 55.335066 seconds (114.95 M allocations: 5.520 GiB, 1.47% gc time, 97.81% compilation time: 45% of which was recompilation)

NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 25-element Vector{Float64}:
0
0
0
0
0
0
0
0
0
0
 ⋮
0
0
0
0
0
0
0
0
0

Solve the PDE

@time sol = NonlinearSolve.solve(prob, NewtonRaphson())

 18.982141 seconds (44.39 M allocations: 1.790 GiB, 0.86% gc time, 70.09% compilation time)

retcode: Stalled
Interpolation: Dict{Symbolics.Num, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}}
ivs: 2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 y
 xdomain:(-0.05:0.1:1.05, -0.05:0.1:1.05)
u: Dict{Symbolics.Num, Matrix{Float64}} with 1 entry:
  u(x, y) => [0.0 -0.00252324 … -0.047774 0.0; -0.00252161 0.00252296 … 0.04751…

Extract data

discrete_x = sol[x]
discrete_y = sol[y]
u_sol = sol[u(x,y)]

12×12 Matrix{Float64}:
  0.0         -0.00252324  -0.00756893  …  -0.0419034  -0.047774   0.0
 -0.00252161   0.00252296   0.00756905      0.0424292   0.0475128  0.0523478
 -0.0075675    0.00756764   0.0227052       0.128131    0.143048   0.158579
 -0.0126098    0.01261      0.0378351       0.213735    0.237969   0.262017
 -0.0176468    0.0176469    0.0529503       0.299786    0.333076   0.364712
 -0.0226742    0.0226742    0.0680376   …   0.387218    0.429839   0.469453
 -0.0276868    0.0276865    0.0830797       0.477544    0.529608   0.534881
 -0.0326795    0.0326789    0.0980596       0.574882    0.676167   0.878422
 -0.0376497    0.0376491    0.112968        0.64799     0.721755   0.730992
 -0.0425993    0.0425989    0.127811        0.725737    0.831872   0.976445
 -0.0475347    0.0475346    0.142609    …   0.78291     0.903549   1.00422
  0.0          0.0524653    0.157391        0.792736    0.995191   0.0

Visualize the solution#

heatmap(
    discrete_x, discrete_y, u_sol,
    xlabel="x values", ylabel="y values", aspect_ratio=:equal,
    title="Steady State Heat Equation", xlims=(0, 1), ylims=(0, 1)
)

Runtime environment#

using Pkg
Pkg.status()

Status `~/work/jl-pde/jl-pde/Project.toml`
  [b0b7db55] ComponentArrays v0.15.30
⌃ [aae7a2af] DiffEqFlux v4.4.0
  [5b8099bc] DomainSets v0.7.16
  [7da242da] Enzyme v0.13.108
  [f6369f11] ForwardDiff v1.3.0
  [d3d80556] LineSearches v7.5.1
⌃ [b2108857] Lux v1.21.1
  [94925ecb] MethodOfLines v0.11.9
  [961ee093] ModelingToolkit v10.31.0
  [315f7962] NeuralPDE v5.20.0
  [8913a72c] NonlinearSolve v4.12.0
⌅ [7f7a1694] Optimization v4.8.0
⌃ [36348300] OptimizationOptimJL v0.4.5
⌃ [42dfb2eb] OptimizationOptimisers v0.3.11
⌃ [500b13db] OptimizationPolyalgorithms v0.3.1
  [1dea7af3] OrdinaryDiffEq v6.103.0
  [91a5bcdd] Plots v1.41.2
  [ce78b400] SimpleUnPack v1.1.0
  [37e2e46d] LinearAlgebra v1.12.0
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

using InteractiveUtils
InteractiveUtils.versioninfo()

Julia Version 1.12.2
Commit ca9b6662be4 (2025-11-20 16:25 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-18.1.7 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 4 default, 1 interactive, 4 GC (on 4 virtual cores)
Environment:
  JULIA_CPU_TARGET = generic;icelake-server,clone_all;znver3,clone_all
  JULIA_CONDAPKG_OFFLINE = true
  JULIA_CONDAPKG_BACKEND = Null
  JULIA_CI = true
  LD_LIBRARY_PATH = /opt/hostedtoolcache/Python/3.13.9/x64/lib
  JULIA_NUM_THREADS = auto

This notebook was generated using Literate.jl.

2D time-independent heat equation

Contents

2D time-independent heat equation#

PDE system#

Visualize the solution#

Runtime environment#