# Problem of convergence for the SH model

**URL:** <https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307>\
**Category:** General Usage\
**Tags:** ode, optimization\
**Created:** [May 17, 2021, 3:05pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307 "2021-05-17T15:05:27Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [May 17, 2021, 3:05pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/1 "2021-05-17T15:05:27Z")

</div>

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:

```julia
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
end
function F_sh(u, p, t)
    @unpack l , L1 = p
	return L1 * u .+ (l .* u)
end
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:

> [CVODES WARNING] CVode  
> Internal t = 0.0446594 and h = 3.61507e-021 are such that t + h = t on the next step. The solver will continue anyway.  
> [CVODES WARNING] CVode  
> The above warning has been issued mxhnil times and will not be issued again for this problem.  
> [CVODES ERROR] CVode  
> 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!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 17, 2021, 3:34pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/2 "2021-05-17T15:34:46Z")

</div>

Try running at a lower tolerance?

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [May 17, 2021, 3:44pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/3 "2021-05-17T15:44:21Z")

</div>

No, how it is done?

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [May 17, 2021, 3:53pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/4 "2021-05-17T15:53:42Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [May 17, 2021, 4:13pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/5 "2021-05-17T16:13:16Z")

</div>

You can try splitODE. See:  
[https://rveltz.github.io/BifurcationKit.jl/dev/tutorialsPD/#Periodic-orbits-from-the-PD-point-(Shooting)-1](https://rveltz.github.io/BifurcationKit.jl/dev/tutorialsPD/#Periodic-orbits-from-the-PD-point-(Shooting)-1)  
for something quite close to your setting

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [May 17, 2021, 4:15pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/6 "2021-05-17T16:15:56Z")

</div>

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 [https://rveltz.github.io/BifurcationKit.jl/dev/tutorials2/#Snaking-in-the-2d-Swift-Hohenberg-equation-1](https://rveltz.github.io/BifurcationKit.jl/dev/tutorials2/#Snaking-in-the-2d-Swift-Hohenberg-equation-1), that’d be

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

```

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

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

```

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [May 17, 2021, 4:17pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/7 "2021-05-17T16:17:02Z")

</div>

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.

Francesco

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [May 17, 2021, 4:32pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/8 "2021-05-17T16:32:06Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [May 17, 2021, 4:33pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/9 "2021-05-17T16:33:47Z")

</div>

there is a sign missing, it should be

```julia
-L1 * u .+ l .* u

```

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [May 17, 2021, 4:41pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/10 "2021-05-17T16:41:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [May 17, 2021, 4:48pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/11 "2021-05-17T16:48:38Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [May 17, 2021, 5:05pm UTC](https://discourse.julialang.org/t/problem-of-convergence-for-the-sh-model/61307/12 "2021-05-17T17:05:58Z")

</div>

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.