Solving the Heat Equation

Author

JMG

using DrWatson
@quickactivate "biomath_des"

using Pkg

Pkg.status()

Status `~/Documents/GitHub/biomath_des/Project.toml`
  [13f3f980] CairoMakie v0.11.4
  [0c46a032] DifferentialEquations v7.12.0
⌅ [5b8099bc] DomainSets v0.6.7
  [634d3b9d] DrWatson v2.13.0
  [61744808] DynamicalSystems v3.2.3
  [94925ecb] MethodOfLines v0.10.4
  [961ee093] ModelingToolkit v8.75.0
  [8913a72c] NonlinearSolve v3.4.0
  [1dea7af3] OrdinaryDiffEq v6.66.0
  [91a5bcdd] Plots v1.39.0
  [f2b01f46] Roots v2.0.22
  [0c5d862f] Symbolics v5.14.0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

using Plots, OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
using Symbolics: scalarize

# Parameters, variables, and derivatives
@parameters t x 
@parameters Dc
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

Differential(x) ∘ Differential(x)

# 1D PDE and boundary conditions
eq  = Dt(u(t, x)) ~ 0.01 * Dxx(u(t, x))
        
bcs = [u(0, x) ~ cos(6*π*x),
        Dx(u(t, 0)) ~ 0.0,
        Dx(u(t, 1)) ~ 0.0]

\[ \begin{align} u\left( 0, x \right) =& \cos\left( 18.85 x \right) \\ \frac{\mathrm{d}}{\mathrm{d}x} u\left( t, 0 \right) =& 0 \\ \frac{\mathrm{d}}{\mathrm{d}x} u\left( t, 1 \right) =& 0 \end{align} \]

# Space and time domains
domains = [t ∈ Interval(0.0, 2.0),
        x ∈ Interval(0.0, 1.0)]

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

# PDE system
@named pdesys = PDESystem(eq, bcs, domains,[t, x],[u(t, x)])

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

# Method of lines discretization
# Need a small dx here for accuracy
dx = 0.01
order = 2
discretization = MOLFiniteDifference([x => dx],t)

MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}(Dict{Num, Float64}(x => 0.01), t, 2, UpwindScheme(1), MethodOfLines.CenterAlignedGrid(), true, false, MethodOfLines.ScalarizedDiscretization(), true, Any[], Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}())

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization);

sol = solve(prob, Tsit5(), saveat=0.05);

sol.u[u(t,x)]

41×101 Matrix{Float64}:
 0.999791     0.982287     0.929776     …  0.982287     0.999791
 0.836977     0.822358     0.778501        0.822358     0.836977
 0.70078      0.688567     0.651931        0.688567     0.70078
 0.586815     0.576578     0.545867        0.576578     0.586815
 0.491369     0.482795     0.457073        0.482795     0.491369
 0.411433     0.404255     0.382722     …  0.404255     0.411433
 0.344492     0.338481     0.320449        0.338481     0.344492
 0.288431     0.283396     0.268293        0.283396     0.288431
 0.241479     0.237263     0.224613        0.237263     0.241479
 0.202156     0.198625     0.188032        0.198625     0.202156
 0.169224     0.166266     0.157394     …  0.166266     0.169224
 0.141643     0.139166     0.131735        0.139166     0.141643
 0.118544     0.116469     0.110246        0.116469     0.118544
 ⋮                                      ⋱               ⋮
 0.00540005   0.00529886   0.00499532      0.00529886   0.00540005
 0.00445675   0.00437213   0.00411827   …  0.00437213   0.00445675
 0.00366777   0.00359695   0.00338451      0.00359695   0.00366777
 0.00300786   0.00294864   0.00277096      0.00294864   0.00300786
 0.00245601   0.00240653   0.00225809      0.00240653   0.00245601
 0.00199476   0.00195339   0.0018293       0.00195339   0.00199476
 0.00160932   0.00157473   0.00147097   …  0.00157473   0.00160932
 0.00128724   0.00125838   0.00117178      0.00125838   0.00128724
 0.00101831   0.000994218  0.000921954     0.000994218  0.00101831
 0.000793874  0.000773753  0.000713393     0.000773753  0.000793874
 0.000606614  0.000589834  0.000539497     0.000589834  0.000606614
 0.000450486  0.000436508  0.000394571  …  0.000436508  0.000450486

discrete_t = sol.t

41-element Vector{Float64}:
 0.0
 0.05
 0.1
 0.15
 0.2
 0.25
 0.3
 0.35
 0.4
 0.45
 0.5
 0.55
 0.6
 ⋮
 1.45
 1.5
 1.55
 1.6
 1.65
 1.7
 1.75
 1.8
 1.85
 1.9
 1.95
 2.0

discrete_x = 0.0:dx:1.0

0.0:0.01:1.0

length(discrete_x)

plt = plot()

for i in eachindex(sol.t)
    plot!(discrete_x, sol.u[u(t,x)][i, :], label="Numerical, t=$(discrete_t[i])"; legend=false)
end
plt

anim = @animate for i in eachindex(discrete_t)
    p1 = plot(discrete_x, sol.u[u(t,x)][i, :], title="u, t=$(discrete_t[i])"; legend=false, xlabel="x",ylabel="u",ylim=(-1,1))
    plot(p1)
end
gif(anim, "plot.gif",fps=100)

[ Info: Saved animation to /Users/jasongraham/Documents/GitHub/biomath_des/notebooks/heat_equation_modeltoolkit/plot.gif