bifurcationkit/BifurcationKit.jl

Issue when implementing the KdV equation

gotencoder opened this issue · 9 comments

Hi there,

I am trying to perform bifurcation analysis on the KdV equation using this package - I am using the pd-1d.jl example script as a template for my implementation. I am having an issue with the function

    newton(probBif, optnewton, normN = norminf)

where an error is thrown ERROR: LoadError: type NamedTuple has no field lx.
I am not sure on the cause of this, I am happy to provide more information/code if that helps. Also to implement something like the KdV is it somewhat straightforward, or are there potential problems to consider?

Many thanks

Hi
I think you forgot to put à parameter axis corresponding to the lens of your problem. Ie there is no getLens(probBof) in your getParams(probBif)

For kdv, it depends 1d or 2d?
Do you want finite difference or FEM which rest on sparse matrices?
You can also use Fourier depending on your BC. In 1d there is a code available by @gibson

Hi there,

Thanks for your reply!

I have added the lx = lx part in the par_br variable and that has seemed to alleviate the problem.
Now the Newton iterations run but the f(x) value gets really large (~5e42) and the code exits with the error 'ERROR: LoadError: SingularException(0)'

I am looking to perform analysis on a 1d kdv (currently I am considering periodic BCs) - any help/resources would be really appreciated! (I am quite new to Julia)

Many thanks!

It should not be...
Can you please show me your problem and parameters?

Absolutely, below is the code I have implemented for a variant of the basic kdv case - I am trying to follow the same structure in the pd-1d.jl script.
All advice is very welcome!

using Revise
using DiffEqOperators, ForwardDiff, DifferentialEquations, Sundials
using BifurcationKit, LinearAlgebra, Plots, SparseArrays, Parameters
const BK = BifurcationKit


norminf(x) = norm(x, Inf)
function Laplacian(N, lx, bc = :Dirichlet)
    hx = 2lx/N
    D2x = CenteredDifference(2, 2, hx, N)
    if bc == :Neumann
        Qx = Neumann0BC(hx)
    elseif bc == :Dirichlet
        Qx = Dirichlet0BC(typeof(hx))
    elseif bc == :Periodic
        Qx = PeriodicBC(typeof(hx))
    end
    D2xsp = sparse(D2x * Qx)[1] |> sparse
end

@views function NL!(dest, u, p, t = 0.)
    N = div(length(u), 2)
    u1 =  u[1:N]
    u2 =  u[N+1:end]

    dest[1:N] .= (1/6) .* (-Laplacian(N, p.lx) * u1 + 9 .* u1 .* Laplacian(N, p.lx) * u1 + Laplacian(N, p.lx, :Neumann) * u1 .* (p.F-1))
    dest[N+1:end] .= (1/6) .* (-Laplacian(N, p.lx) * u2 + 9 .* u1 .* Laplacian(N, p.lx) * u2 + Laplacian(N, p.lx, :Neumann) * u2 .* (p.F-1))


    return dest
end

function Fbr!(f, u, p, t = 0.)
    NL!(f, u, p)
    mul!(f, p.Δ, u,1,1)
    f
end

NL(u, p) = NL!(similar(u), u, p)
Fbr(x, p, t = 0.) = Fbr!(similar(x), x, p)
dNL(x, p, dx) = ForwardDiff.derivative(t -> NL(x .+ t .* dx, p), 0.)

function dFbr!(f, x, p, dx)
    mul!(f, p.Δ, dx)

    nl = dNL(x, p, dx)
    f .= f .+ nl
end
dFbr(x, p, dx) = dFbr!(similar(dx), x, p, dx)

Jbr(x, p) = sparse(ForwardDiff.jacobian(x -> Fbr(x, p), x))

#Set params

N = 100
n = 2N
lx = 3pi /2
X = LinRange(-lx,lx, N)
Δ = Laplacian(N, lx, :Neumann)
par_br = (lx = lx, F = 1.1, Δ = blockdiag(Δ, Δ))

u0 = cos.(2X)
solc0 = vcat(u0, u0)


probBif = BifurcationProblem(Fbr, solc0, par_br, (@lens _.F) ;J = Jbr,
		recordFromSolution = (x, p) -> norm(x, Inf),
		plotSolution = (x, p; kwargs...) -> plot!(x[1:end÷2];label="",ylabel ="u", kwargs...))


#So far the code runs without an error
############################################

eigls = EigArpack(0.5, :LM)
optnewton = NewtonPar(eigsolver = eigls, verbose=true, maxIter = 3200, tol=1e-9)

out = @time newton(probBif, optnewton, normN = norminf)

The code looks fine for BifurcationKit. I guess your implementation of kdv is not good, your initial guess is not good enough, etc.

Should we close this?

Thanks for the guidance, sure happy to close this.
Is there an email/method to contact you rather than opening up a new 'issue'?

Thanks again

romain veltz inria fr