2D 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 Plots
[ Info: Precompiling IJuliaExt [2f4121a4-3b3a-5ce6-9c5e-1f2673ce168a]
Setup variables and differential operators
@variables t x y
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Differential(y) ∘ Differential(y)
PDE equation
eq = Dt(u(t, x, y)) ~ Dxx(u(t, x, y)) + Dyy(u(t, x, y))
\[ \begin{equation}
\frac{\mathrm{d}}{\mathrm{d}t} u\left( t, x, y \right) = \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} u\left( t, x, y \right) + \frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} u\left( t, x, y \right)
\end{equation}
\]
Boundary conditions
bcs = [
u(0, x, y) ~ 0,
u(t, 0, y) ~ x * y,
u(t, 1, y) ~ x * y,
u(t, x, 0) ~ x * y,
u(t, x, 1) ~ x * y,
]
\[\begin{split} \begin{align}
u\left( 0, x, y \right) =& 0 \\
u\left( t, 0, y \right) =& x y \\
u\left( t, 1, y \right) =& x y \\
u\left( t, x, 0 \right) =& x y \\
u\left( t, x, 1 \right) =& x y
\end{align}
\end{split}\]
Space and time domains
domains = [
t ∈ Interval(0.0, 1.0),
x ∈ Interval(0.0, 1.0),
y ∈ Interval(0.0, 1.0)
]
3-element Vector{Symbolics.VarDomainPairing}:
Symbolics.VarDomainPairing(t, 0.0 .. 1.0)
Symbolics.VarDomainPairing(x, 0.0 .. 1.0)
Symbolics.VarDomainPairing(y, 0.0 .. 1.0)
PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])
\[ \begin{align}
\frac{\mathrm{d}}{\mathrm{d}t} u\left( t, x, y \right) =& \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} u\left( t, x, y \right) + \frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} u\left( t, x, y \right)
\end{align}
\]
discretize the PDE system
N = 20
discretization = MOLFiniteDifference([x=>N, y=>N], t, approx_order=2, grid_align=MethodOfLines.EdgeAlignedGrid())
prob = discretize(pdesys, discretization)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 400-element Vector{Float64}:
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
⋮
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
-0.0
Solve the PDE
sol = solve(prob, KenCarp4(), saveat=0.01)
retcode: Success
Interpolation: Dict{Num, Interpolations.GriddedInterpolation{Float64, 3, Array{Float64, 3}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}}}
t: 101-element Vector{Float64}:
0.0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0.11
0.12
⋮
0.89
0.9
0.91
0.92
0.93
0.94
0.95
0.96
0.97
0.98
0.99
1.0ivs: 3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
t
x
ydomain:([0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09 … 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], -0.025:0.05:1.025, -0.025:0.05:1.025)
u: Dict{Num, Array{Float64, 3}} with 1 entry:
u(t, x, y) => [0.0 0.0 … 0.0 0.0; 0.0 1.6849e-10 … -0.0214844 0.0; … ; 0.0 -0…
Extract data
discrete_x = sol[x]
discrete_y = sol[y]
discrete_t = sol[t]
solu = sol[u(t, x, y)]
101×22×22 Array{Float64, 3}:
[:, :, 1] =
0.0 0.0 0.0 … 0.0 0.0 0.0
0.0 1.6849e-10 5.24652e-10 -0.0148733 -0.0214844 0.0
0.0 -9.69678e-8 -3.62435e-7 -0.0176809 -0.0225107 0.0
0.0 -2.79302e-6 -9.19576e-6 -0.0189859 -0.0229673 0.0
0.0 -1.59429e-5 -5.04371e-5 -0.0198141 -0.023254 0.0
0.0 -4.46441e-5 -0.000138676 … -0.0204107 -0.0234577 0.0
0.0 -8.88607e-5 -0.00027295 -0.0208908 -0.0236212 0.0
0.0 -0.000142968 -0.0004361 -0.0212853 -0.023755 0.0
0.0 -0.000201733 -0.000612571 -0.0216081 -0.0238641 0.0
0.0 -0.000260136 -0.000787453 -0.0218747 -0.0239541 0.0
0.0 -0.000313856 -0.000947999 … -0.0220977 -0.0240293 0.0
0.0 -0.000361957 -0.00109159 -0.0222817 -0.0240914 0.0
0.0 -0.000404475 -0.00121839 -0.0224324 -0.0241421 0.0
⋮ ⋱ ⋮
0.0 -0.000624994 -0.00187498 -0.023125 -0.024375 0.0
0.0 -0.000624991 -0.00187497 … -0.023125 -0.024375 0.0
0.0 -0.000624989 -0.00187497 -0.023125 -0.024375 0.0
0.0 -0.000624987 -0.00187496 -0.023125 -0.024375 0.0
0.0 -0.000624986 -0.00187496 -0.023125 -0.024375 0.0
0.0 -0.000624986 -0.00187496 -0.023125 -0.024375 0.0
0.0 -0.000624987 -0.00187496 … -0.023125 -0.024375 0.0
0.0 -0.000624988 -0.00187496 -0.023125 -0.024375 0.0
0.0 -0.00062499 -0.00187497 -0.023125 -0.024375 0.0
0.0 -0.000624991 -0.00187497 -0.023125 -0.024375 0.0
0.0 -0.000624993 -0.00187498 -0.023125 -0.024375 0.0
0.0 -0.000624994 -0.00187498 … -0.023125 -0.024375 0.0
[:, :, 2] =
0.0 -0.0 … -0.0 -0.0 0.08
1.6849e-10 -1.94823e-10 0.0148733 0.0214844 0.0285365
-9.69678e-8 8.90186e-8 0.0176808 0.0225107 0.027491
-2.79302e-6 2.75146e-6 0.0189857 0.0229672 0.0270343
-1.59429e-5 1.58755e-5 0.0198137 0.0232539 0.0267459
-4.46441e-5 4.45774e-5 … 0.0204104 0.0234576 0.0265428
-8.88607e-5 8.88119e-5 0.0208906 0.0236211 0.026379
-0.000142968 0.000142939 0.0212852 0.023755 0.026245
-0.000201733 0.000201721 0.021608 0.0238641 0.0261359
-0.000260136 0.000260135 0.0218747 0.0239541 0.026046
-0.000313856 0.000313861 … 0.0220977 0.0240293 0.0259707
-0.000361957 0.000361966 0.0222818 0.0240914 0.0259086
-0.000404475 0.000404484 0.0224324 0.0241421 0.0258579
⋮ ⋱ ⋮
-0.000624994 0.000624994 0.023125 0.024375 0.025625
-0.000624991 0.000624991 … 0.023125 0.024375 0.025625
-0.000624989 0.000624989 0.023125 0.024375 0.025625
-0.000624987 0.000624987 0.023125 0.024375 0.025625
-0.000624986 0.000624986 0.023125 0.024375 0.025625
-0.000624986 0.000624986 0.023125 0.024375 0.025625
-0.000624987 0.000624987 … 0.023125 0.024375 0.025625
-0.000624988 0.000624988 0.023125 0.024375 0.025625
-0.00062499 0.00062499 0.023125 0.024375 0.025625
-0.000624991 0.000624991 0.023125 0.024375 0.025625
-0.000624993 0.000624993 0.023125 0.024375 0.025625
-0.000624994 0.000624994 … 0.023125 0.024375 0.025625
[:, :, 3] =
0.0 -0.0 … -0.0 -0.0 0.24
5.24652e-10 -6.03652e-10 0.0446199 0.0644531 0.0856094
-3.62435e-7 3.38588e-7 0.053043 0.0675322 0.0824728
-9.19576e-6 9.07112e-6 0.0569624 0.0689035 0.0811012
-5.04371e-5 5.02353e-5 0.0594553 0.0697665 0.080233
-0.000138676 0.000138476 … 0.0612535 0.0703802 0.0796208
-0.00027295 0.000272804 0.0626984 0.0708723 0.0791279
-0.0004361 0.000436013 0.0638833 0.0712743 0.0787255
-0.000612571 0.000612535 0.0648506 0.0716012 0.0783988
-0.000787453 0.000787452 0.0656481 0.0718703 0.0781299
-0.000947999 0.000948015 … 0.0663142 0.0720951 0.077905
-0.00109159 0.00109162 0.0668636 0.0722803 0.0777197
-0.00121839 0.00121842 0.067313 0.0724316 0.0775684
⋮ ⋱ ⋮
-0.00187498 0.00187498 0.0693749 0.073125 0.076875
-0.00187497 0.00187497 … 0.0693749 0.073125 0.076875
-0.00187497 0.00187497 0.0693749 0.073125 0.076875
-0.00187496 0.00187496 0.0693749 0.073125 0.076875
-0.00187496 0.00187496 0.0693749 0.073125 0.076875
-0.00187496 0.00187496 0.0693749 0.073125 0.076875
-0.00187496 0.00187496 … 0.0693749 0.073125 0.076875
-0.00187496 0.00187496 0.0693749 0.073125 0.076875
-0.00187497 0.00187497 0.0693749 0.073125 0.076875
-0.00187497 0.00187497 0.0693749 0.073125 0.076875
-0.00187498 0.00187498 0.0693749 0.073125 0.076875
-0.00187498 0.00187498 … 0.0693749 0.073125 0.076875
;;; …
[:, :, 20] =
0.0 -0.0 -0.0 … -0.0 -0.0 2.96
-0.0148733 0.0148733 0.0446199 0.746433 0.863582 0.986738
-0.0176809 0.0176808 0.053043 0.808188 0.885656 0.964307
-0.0189859 0.0189857 0.0569624 0.828028 0.892483 0.957534
-0.0198141 0.0198137 0.0594553 0.837704 0.895821 0.954154
-0.0204107 0.0204104 0.0612535 … 0.842912 0.897578 0.952423
-0.0208908 0.0208906 0.0626984 0.846325 0.89874 0.951256
-0.0212853 0.0212852 0.0638833 0.848654 0.899529 0.950465
-0.0216081 0.021608 0.0648506 0.850204 0.90005 0.949948
-0.0218747 0.0218747 0.0656481 0.851317 0.900424 0.949576
-0.0220977 0.0220977 0.0663142 … 0.852228 0.900732 0.949267
-0.0222817 0.0222818 0.0668636 0.852927 0.900968 0.949031
-0.0224324 0.0224324 0.067313 0.853456 0.901146 0.948853
⋮ ⋱ ⋮
-0.023125 0.023125 0.0693749 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 … 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 … 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 0.855625 0.901875 0.948125
-0.023125 0.023125 0.0693749 … 0.855625 0.901875 0.948125
[:, :, 21] =
0.0 -0.0 -0.0 … -0.0 -0.0 3.12
-0.0214844 0.0214844 0.0644531 0.863582 0.937203 1.01289
-0.0225107 0.0225107 0.0675322 0.885656 0.945115 1.0048
-0.0229673 0.0229672 0.0689035 0.892483 0.947428 1.00258
-0.023254 0.0232539 0.0697665 0.895821 0.94859 1.00138
-0.0234577 0.0234576 0.0703802 … 0.897578 0.949174 1.00083
-0.0236212 0.0236211 0.0708723 0.89874 0.949569 1.00043
-0.023755 0.023755 0.0712743 0.899529 0.949837 1.00016
-0.0238641 0.0238641 0.0716012 0.90005 0.950012 0.999987
-0.0239541 0.0239541 0.0718703 0.900424 0.950136 0.999864
-0.0240293 0.0240293 0.0720951 … 0.900732 0.950241 0.999759
-0.0240914 0.0240914 0.0722803 0.900968 0.95032 0.999679
-0.0241421 0.0241421 0.0724316 0.901146 0.95038 0.999619
⋮ ⋱ ⋮
-0.024375 0.024375 0.073125 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 … 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 … 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 0.901875 0.950625 0.999375
-0.024375 0.024375 0.073125 … 0.901875 0.950625 0.999375
[:, :, 22] =
0.0 0.08 0.24 0.4 … 2.8 2.96 3.12 0.0
0.0 0.0285365 0.0856094 0.142682 0.95503 0.986738 1.01289 0.0
0.0 0.027491 0.0824728 0.137454 0.922459 0.964307 1.0048 0.0
0.0 0.0270343 0.0811012 0.135162 0.911921 0.957534 1.00258 0.0
0.0 0.0267459 0.080233 0.133705 0.906642 0.954154 1.00138 0.0
0.0 0.0265428 0.0796208 0.132676 … 0.90384 0.952423 1.00083 0.0
0.0 0.026379 0.0791279 0.131849 0.901973 0.951256 1.00043 0.0
0.0 0.026245 0.0787255 0.131178 0.900698 0.950465 1.00016 0.0
0.0 0.0261359 0.0783988 0.130635 0.899853 0.949948 0.999987 0.0
0.0 0.026046 0.0781299 0.13019 0.899245 0.949576 0.999864 0.0
0.0 0.0259707 0.077905 0.129818 … 0.898744 0.949267 0.999759 0.0
0.0 0.0259086 0.0777197 0.129512 0.898359 0.949031 0.999679 0.0
0.0 0.0258579 0.0775684 0.129263 0.898069 0.948853 0.999619 0.0
⋮ ⋱ ⋮
0.0 0.025625 0.076875 0.128125 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 … 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 … 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 0.896875 0.948125 0.999375 0.0
0.0 0.025625 0.076875 0.128125 … 0.896875 0.948125 0.999375 0.0
Animate the solution
anim = @animate for k in eachindex(discrete_t)
heatmap(solu[k, 2:end-1, 2:end-1], title="u @ t=$(discrete_t[k])", aspect_ratio=:equal)
end
mp4(anim,"2Dheat.mp4", fps = 8)
[ Info: Saved animation to /tmp/docs/pde/2Dheat.mp4
Runtime environment#
using Pkg
Pkg.status()
Status `/tmp/docs/pde/Project.toml`
⌅ [5b8099bc] DomainSets v0.6.7
[94925ecb] MethodOfLines v0.11.0
[961ee093] ModelingToolkit v9.15.0
[8913a72c] NonlinearSolve v3.12.4
[1dea7af3] OrdinaryDiffEq v6.80.1
[91a5bcdd] Plots v1.40.4
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`
using InteractiveUtils
InteractiveUtils.versioninfo()
Julia Version 1.10.3
Commit 0b4590a5507 (2024-04-30 10:59 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
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 4 default, 0 interactive, 2 GC (on 4 virtual cores)
Environment:
JULIA_CI = true
JULIA_NUM_THREADS = auto
JULIA_CONDAPKG_BACKEND = Null
JULIA_PATH = /usr/local/julia/
JULIA_DEPOT_PATH = /srv/juliapkg/