# Stability of steady-state heat equation example for MethodOfLines.jl
**URL:** https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651
**Category:** Modelling & Simulations
**Tags:** pde
**Created:** [May 10, 2023, 11:25pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651 "2023-05-10T23:25:44Z")
**Posts on this page:** 17
**Page:** 1
### Author:  [@thomasteisberg](https://discourse.julialang.org/u/thomasteisberg)
#### Post date: [May 10, 2023, 11:25pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/1 "2023-05-10T23:25:45Z")
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](https://docs.sciml.ai/MethodOfLines/stable/tutorials/heatss/)
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`):
```julia
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.
---
### Author:  [@xtalax](https://discourse.julialang.org/u/xtalax)
#### Post date: [May 11, 2023, 1:46pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/2 "2023-05-11T13:46:51Z")
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](https://docs.sciml.ai/MethodOfLines/stable/nonuniform/)?
---
### Author:  [@thomasteisberg](https://discourse.julialang.org/u/thomasteisberg)
#### Post date: [May 11, 2023, 4:34pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/3 "2023-05-11T16:34:12Z")
Hmm still behaves weirdly past 12 points in this example.
```julia
cheb_n = 14
discx = chebyspace(cheb_n, domains[1])
discy = chebyspace(cheb_n, domains[2])
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!
---
### Author:  [@xtalax](https://discourse.julialang.org/u/xtalax)
#### Post date: [May 11, 2023, 4:36pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/4 "2023-05-11T16:36:06Z")
Can you link this thread in an issue on the repo please?
---
### Author:  [@thomasteisberg](https://discourse.julialang.org/u/thomasteisberg)
#### Post date: [May 11, 2023, 5:17pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/5 "2023-05-11T17:17:01Z")
Sure. Here it is: [Unstable steady-state solutions for smaller grid sizes in heat equation example · Issue #274 · SciML/MethodOfLines.jl · GitHub](https://github.com/SciML/MethodOfLines.jl/issues/274)
---
### Author:  [@parb](https://discourse.julialang.org/u/parb)
#### Post date: [October 29, 2023, 9:32am UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/6 "2023-10-29T09:32:24Z")
A solution here would be really great! I’m working on a problem now where I would’ve loved to have used this MethodOfLines.jl. (I use the Laplacian on a cylindrical geometry, and with slightly different boundary conditions, but have the exact same problem occur for dx = dy \< 0.1).
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [August 31, 2024, 4:03pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/8 "2024-08-31T16:03:46Z")
Just stumbled upon the problem ☹
Step of 0.1 works fine, step 0.05 diverges and stops on MaxIters, and step 0.025 or less appears to hang Julia, which doesn’t then react to Ctrl-C.
Are there any workarounds known (using some other algorithm instead of `NewtonRaphson` ??).
Alternative question: If I have a relatively simple 2D stationary diffusion problem, what is currently the way to go?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 31, 2024, 4:05pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/9 "2024-08-31T16:05:21Z")
> [@Eben60](#):
>
> Are there any workarounds known (using some other algorithm instead of `NewtonRaphson` ??).
Use RobustMultiNewton()?
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [August 31, 2024, 4:06pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/10 "2024-08-31T16:06:17Z")
> [@ChrisRackauckas](#):
>
> RobustMultiNewton()
were you specially sitting here waiting for my question? 😀
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 31, 2024, 4:06pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/11 "2024-08-31T16:06:41Z")
No, just happened to be here… all the time.
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [August 31, 2024, 4:14pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/12 "2024-08-31T16:14:25Z")
it’s the exactly same with `RobustMultiNewton` : Diverges at step 0.05 (exactly the same picture as here [MethodOfLines - numerical instability for example in docs](https://discourse.julialang.org/t/methodoflines-numerical-instability-for-example-in-docs/106359)), and hangs at 0.025 (Julia continues to run in the background even is I “trash” it in the VSCode)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 31, 2024, 4:26pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/13 "2024-08-31T16:26:48Z")
That might need an issue
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [August 31, 2024, 4:31pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/14 "2024-08-31T16:31:04Z")
The issue was already linked above - [Unstable steady-state solutions for smaller grid sizes in heat equation example · Issue #274 · SciML/MethodOfLines.jl · GitHub](https://github.com/SciML/MethodOfLines.jl/issues/274)
I’m adding link to my posts there
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 31, 2024, 4:40pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/15 "2024-08-31T16:40:17Z")
Thanks
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [August 31, 2024, 9:47pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/16 "2024-08-31T21:47:33Z")
As the problem in the example is actually linear, I wondered why `NonlinearSolve` is to be used, and tried `LinearSolve.solve(prob, nothing)`. That worked for the step of 0.1 and diverged for 0.05 just the same.
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [September 1, 2024, 6:27am UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/17 "2024-09-01T06:27:29Z")
That probably means that it’s not the solver’s fault right? (Because 3 different solvers experience the same issue) So perhaps something with the discretization goes awry?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 1, 2024, 2:59pm UTC](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/18 "2024-09-01T14:59:31Z")
Yes. Or this could be a property of the discretization on this equation. Whether stepping is possible to be stable is a property of using the right spatial semi-discretization, and I haven’t looked into that detail here.