I am using the steady state heat equation example in the MethodOfLines.jl package.
When I reduce the step side below 0.1 I have a numerical blow-up of the solution.
I do not know if this is only on my end or other users have also faced this error.
To me I do not see a reason for numerical instability of this case. I have attached my code below:
using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve
@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ 0
bcs = [u(0, y) ~ x * y,
u(1, y) ~ x * y,
u(x, 0) ~ x * y,
u(x, 1) ~ x * y]
# Space and time domains
domains = [x ∈ Interval(0.0, 1.0),
y ∈ Interval(0.0, 1.0)]
@named pdesys = PDESystem([eq], bcs, domains, [x, y], [u(x, y)])
dx = 0.05
dy = 0.05
# Note that we pass in `nothing` for the time variable `t` here since we
# are creating a stationary problem without a dependence on time, only space.
discretization = MOLFiniteDifference([x => dx, y => dy], nothing, approx_order=2)
prob = discretize(pdesys, discretization)
sol = NonlinearSolve.solve(prob, NewtonRaphson())
u_sol = sol[u(x, y)]
using Plots
heatmap(sol[x], sol[y], u_sol, xlabel="x values", ylabel="y values",
title="Steady State Heat Equation")
This is the solution I get:
