# Stability of steady-state heat equation example for MethodOfLines.jl

The steady-state heat equation example from the MethodOfLines.jl documentation seems to break for smaller step sizes. The example I’m referring to is here: Steady State Heat Equation - No Time Dependence - NonlinearProblem · MethodOfLines.jl

Output for `dx == dy == 0.1` (as in the original tutorial):

Output for `dx == dy == 0.09`:

Output for `dx == dy == 0.05` (all `NaN`s on the interior):

The return code in every `PDESolution` is `MaxIters` (including the `0.1` case), so I also tried significantly increasing `maxiters` from the default. This doesn’t seem to do much:

Is there some reason why this is expected behavior?

Minimal example (essentially identical to the example in the docs, except that I’m using Makie for plotting and varying `dx`/`dy` and `maxiters`):

``````using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve
using CairoMakie

@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, dy = 0.1,  0.1   # works great
#dx, dy = 0.09, 0.09 # instability near (1, 1)
#dx, dy = 0.08, 0.08 # unstable region grows
#dx, dy = 0.05, 0.05 # all NaN on the interior of the domain

iters_power = 5
maxiters = Int(10^iters_power)

# 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(), maxiters=maxiters)

begin
fig = Figure(resolution=(1000, 800))
ax = Axis(fig[1, 1], title="dx = \$dx, dy = \$dy, maxiters = 1e\$iters_power")
h = heatmap!(ax, sol[x], sol[y], sol[u(x, y)], colorrange=(0, 1))
cb = Colorbar(fig[1, 2], h, label="u(x, y)")
fig
end

save("heat_ss_\$(dx)_1e\$(iters_power).png", fig)
``````

Thanks in advance for taking a look!

Edit: Updated with more examples.

1 Like

Looks like a problem with nonuniform grids and nonlinear problem that has gone undetected so far, I will do some more benchmarking to determine the cause.

Can you see if this is present using chebyspace?

1 Like

Hmm still behaves weirdly past 12 points in this example.

``````cheb_n = 14
discx = chebyspace(cheb_n, domains)
discy = chebyspace(cheb_n, domains)
discretization = MOLFiniteDifference([discx, discy], nothing, approx_order=2)

prob = discretize(pdesys, discretization)
sol = NonlinearSolve.solve(prob, NewtonRaphson(), maxiters=maxiters)

begin
fig = Figure(resolution=(1000, 800))
ax = Axis(fig[1, 1], title="chebyspace points per axis = \$cheb_n, maxiters = 1e\$iters_power")
h = heatmap!(ax, sol[x], sol[y], sol[u(x, y)], colorrange=(0, 1))
cb = Colorbar(fig[1, 2], h, label="u(x, y)")
fig
end

save("heat_ss_cheb\$(cheb_n)_1e\$(iters_power).png", fig)
``````

Thanks for taking a look!

1 Like