swh:1:snp:96df3cd8f41a04650ca4e93b2610f839c02d2899
Tip revision: 60df0a6859f42cad7170a8560c768ff239190d3a authored by Software Heritage on 19 March 2019, 08:28:07 UTC
hal: Deposit 241 in collection hal
hal: Deposit 241 in collection hal
Tip revision: 60df0a6
Langmuir.jl
using Revise
using PseudoArcLengthContinuation, LinearAlgebra, Plots, ApproxFun
const Cont = PseudoArcLengthContinuation
####################################################################################################
# specific methods for ApproxFun
import Base: length, eltype, copyto!
import LinearAlgebra: norm, dot
eltype(x::ApproxFun.Fun) = eltype(x.coefficients)
length(x::ApproxFun.Fun) = length(x.coefficients)
norm(x::ApproxFun.Fun, p::Real) = (@show p;norm(x.coefficients, p))
norm(x::Array{Fun, 1}, p::Real) = (@show p;norm(x[3].coefficients, p))
norm(x::Array{Fun{Chebyshev{Segment{Float64}, Float64}, Float64, Array{Float64, 1}}, 1}, p::Real) = (@show p;norm(x[3].coefficients, p))
dot(x::ApproxFun.Fun, y::ApproxFun.Fun) = sum(x * y)
dot(x::Array{Fun{Chebyshev{Segment{Float64}, Float64}, Float64, Array{Float64, 1}}, 1}, y::Array{Fun{Chebyshev{Segment{Float64}, Float64}, Float64, Array{Float64, 1}}, 1}) = sum(x[3]*y[3])
copyto!(x::ApproxFun.Fun, y::ApproxFun.Fun) = (x.coefficients = y.coefficients)
####################################################################################################
function F(c, V, μ = 0.5, c0 = -0.9)
# plot(c) |> display
return [Fun(c(0), domain(c)) - c0,
Fun((∂ * c)(L), domain(c)),
Fun((Δ * c)(0), domain(c)),
Fun((Δ * c)(L), domain(c)),
-(c'' - c^3 + c - μ * ζ )'' - V * c']
end
function Jac(c, V, μ = 0.5, c0 = -0.9)
return [Evaluation(c.space, 0.),
Evaluation(c.space, L, 1),
Evaluation(c.space, 0, 2),
Evaluation(c.space, L, 2),
-Δ * (Δ - 3c^2 + I) - V * ∂]
end
const L = 100.
sol = Fun(x -> -1.0 , Interval(0.0, L))
xs = 10.
ls = 2.
const ζ = -1/2Fun(x -> 1 + tanh((x-xs)/ls), Interval(0.0, L))
const Δ = Derivative(sol.space, 2)
const ∂ = Derivative(sol.space, 1)
F(sol,0.07) |> norm
opt_new = Cont.NewtonPar(tol = 1e-12, verbose = true)
out, hist, flag = @time Cont.newton(
x -> F(x, 0.07, 0., -0.9),
x -> Jac(x, 0.07, 0., -0.9),
sol, opt_new, normN = x -> norm(x, Inf64))
plot(sol)
opts_br0 = ContinuationPar(dsmin = 0.0001, dsmax = 0.01, ds= 0.001, a = 1.0, pMax = 0.5, theta = 0.5, secant = true, plot_every_n_steps = 3, newtonOptions = NewtonPar(tol = 1e-6, maxIter = 50, verbose = true), doArcLengthScaling = false)
opts_br0.newtonOptions.damped = false
opts_br0.detect_fold = true
opts_br0.maxSteps = 14300
br, u1 = @time Cont.continuation(
(x, p) -> F(x, 0.07, p, -0.9),
(x, p) -> Jac(x, 0.07, p, -0.9),
out, 0.0, opts_br0,
plot = true,
# finaliseSolution = finalise_solution,
plotsolution = (x;kwargs...) -> plot!(x, subplot = 4, label = "l = $(length(x))"))