Problem of convergence for the SH model

Hi everyone,

I’m trying to solving the Swift Hohenberg 2D model as a ODE problem, but it works only on a short time range and reducing the resolution of my grid.
I believe that the problem is not only an optimization problem, but what do you think about? Here’s my code:

using DifferentialEquations
using DiffEqOperators, Setfield, Parameters
using BifurcationKit, LinearAlgebra, Plots, SparseArrays
const BK = BifurcationKit

Nx = 128
Ny = 128
lx = 4pi
ly = 4pi

# we use DiffEqOperators to compute the Laplacian operator
function Laplacian2D(Nx, Ny, lx, ly)
	hx = 2*lx/Nx
	hy = 2*ly/Ny
	D2x = CenteredDifference(2, 2, hx, Nx)
	D2y = CenteredDifference(2, 2, hy, Ny)
	Qx = Neumann0BC(hx)
	Qy = Neumann0BC(hy)
	A = kron(sparse(I, Ny, Ny), sparse(D2x * Qx)[1]) + kron(sparse(D2y * Qy)[1], sparse(I, Nx, Nx))
	return A
function F_sh(u, p, t)
    @unpack l , L1 = p
	return L1 * u .+ (l .* u)
X = -lx .+ 2lx/(Nx) * collect(0:Nx-1)
Y = -ly .+ 2ly/(Ny) * collect(0:Ny-1)

# define parameters for the PDE
Δ = Laplacian2D(Nx, Ny, lx, ly)

L1 = (I + Δ)^2;

#l= -0.1

p = (l= -0.1, L1=L1)
u0=   rand(Nx,Ny) |> vec
tspan = (0.0,100.0)
prob = ODEProblem(F_sh,u0,tspan,p)

 sol =solve(prob, alg_hints=[:stiff],save_everystep=false, save_start=false)

And the error is returned:

Internal t = 0.0446594 and h = 3.61507e-021 are such that t + h = t on the next step. The solver will continue anyway.
The above warning has been issued mxhnil times and will not be issued again for this problem.
At t = 0.0446594 and h = 3.61829e-023, the corrector convergence test failed repeatedly or with |h| = hmin.

Thank you so much for the help, every comment is welcome!

Try running at a lower tolerance?

No, how it is done?

I tried now changing the reltol and abstol values, but the error is still returned… maybe the problem is in the definition of the function?

You can try splitODE. See:
for something quite close to your setting

It looks like you’re solving just the linear terms of the Swift-Hohenberg equation and have left out the nonlinear damping term (e.g. -u^3). I’d guess the linear equation is unstable for the parameters you’ve chosen. Could it be that the solution is blowing up, causing the time step to shrink in order to meet an absolute tolerance condition?

First thing I’d do is put back in the nonlinear terms and see if the solution converges to something sensible. From Snaking in the 2d Swift-Hohenberg equation · Bifurcation Analysis in Julia, that’d be

function F_sh(u, p)
	@unpack l, ν, L1 = p
	return -L1 * u .+ (l .* u .+ ν .* u.^2 .- u.^3)

And maybe while you’re at it change the domain size and discretization to match theirs, just to rule out that difference

Nx = 151
Ny = 100
lx = 4*2pi
ly = 2*2pi/sqrt(3)

I saw that if I use only L1= Δ the convergence is respected. So the problem begin when I write L1=Δ^2.
There is another way to write the operator L1 successfully?

Thank you very much.


Yes, for the moment I’m interested to study SH with only the linear term. Here only the solution u=0 is stable, then starting from random initial condition I expect to find the system at u=0. But maybe I’m doing something wrong?

there is a sign missing, it should be

-L1 * u .+ l .* u

Nice spot. That would explain blowing up on such a short time scale.

I’m so sorry for the very stupid mistake.
Thank you so much for the help, probably I’d have never seen it.

1 Like

Sometimes it takes someone else’s eyes! That seems to be the nature of software. And personally, I learn things from reading & participating in these discussions, even if the final answer is a typo.

1 Like