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