2D time-independent heat equation

From the tutorial

• https://docs.sciml.ai/MethodOfLines/stable/tutorials/heatss/

• https://docs.sciml.ai/MethodOfLines/stable/tutorials/heat/

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

@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

prob = let dx=0.1
    discretization = MOLFiniteDifference([x=>dx, y=>dx], nothing, approx_order=2)
    prob = discretize(pdesys, discretization)
end

NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 9-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

Solve the PDE

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

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.0:0.1:1.0, 0.0:0.1:1.0)
u: Dict{Symbolics.Num, Matrix{Float64}} with 1 entry:
  u(x, y) => [0.0 0.0 … 0.0 0.0; 0.0 0.01 … 0.09 0.1; … ; 0.0 0.09 … 0.809998 0…

Extract data

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

11×11 Matrix{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.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09      0.1
 0.0  0.02  0.04  0.06  0.08  0.1   0.12  0.14  0.16  0.18      0.2
 0.0  0.03  0.06  0.09  0.12  0.15  0.18  0.21  0.24  0.27      0.3
 0.0  0.04  0.08  0.12  0.16  0.2   0.24  0.28  0.32  0.36      0.4
 0.0  0.05  0.1   0.15  0.2   0.25  0.3   0.35  0.4   0.45      0.5
 0.0  0.06  0.12  0.18  0.24  0.3   0.36  0.42  0.48  0.540001  0.6
 0.0  0.07  0.14  0.21  0.28  0.35  0.42  0.49  0.56  0.629998  0.7
 0.0  0.08  0.16  0.24  0.32  0.4   0.48  0.56  0.64  0.720002  0.8
 0.0  0.09  0.18  0.27  0.36  0.45  0.54  0.63  0.72  0.809998  0.9
 0.0  0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9       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)
)

---

This notebook was generated using Literate.jl.