swh:1:snp:96df3cd8f41a04650ca4e93b2610f839c02d2899
Raw File
Tip revision: 60df0a6859f42cad7170a8560c768ff239190d3a authored by Software Heritage on 19 March 2019, 08:28:07 UTC
hal: Deposit 241 in collection hal
Tip revision: 60df0a6
SHpde_snaking.jl
using Revise
using SparseArrays, LinearAlgebra, DiffEqOperators
using PseudoArcLengthContinuation
using Plots
const Cont = PseudoArcLengthContinuation
################################################################################
# case of the SH equation

nx = 400; Lx = 50; hx = 2*Lx/nx
const X = -Lx .+ 2Lx/nx*(1:nx) |> collect
BC = :periodic
Dxx = sparse(DerivativeOperator{Float64}(2, 2, hx, nx, BC, BC))
Lsh = -(I+Dxx)^2

function R_SH(u, p, b, L1)
	out = similar(u)
	out .= L1 * u .- p .* u .+ b .* u.^3 - u.^5
end

SHJac = (u, v, p, b, L1) -> L1*v .+ (-p .+ 3*b .* u.^2 .- 5*u.^4) .* v
Jac_sp = (u, p, b, L1) -> L1 + spdiagm(0 => -p .+ 3*b .* u.^2 .- 5*u.^4)

Fpde = (u, p)-> R_SH(u, p, 2., Lsh)
# jacobian version with full matrix, waste of ressources!!
Jac_fd(u0, alpha) = Cont.finiteDifferences(u->Fpde(u, alpha), u0)

sol0 = 1.65cos.(X) .* exp.(-X.^2/(2*5^2))
	opt_new = Cont.NewtonPar(verbose = true, tol = 1e-11)
	# allocations 26.47k, 0.038s, tol = 1e-10
	sol1, hist, flag = @time Cont.newton(
				u ->   R_SH(u, 0.7, 2., Lsh),
				u -> Jac_sp(u, 0.7, 2., Lsh),
				sol0, opt_new)
	Plots.plot(X,sol1)

if 1==0
	# this is for testing
	# use of Sparse Jacobian
	sol_, hist, flag = @time Cont.newton(
				x->  R_SH(x, 0.7, 2., Lsh),
				u->Jac_sp(u, 0.7, 2., Lsh),
				u0, opt_new)

	# use of FD Jacobian
	sol_, hist, flag = @time Cont.newton(
				x->R_SH(x, 0.7, 2., Lsh),
				u->Jac_fd(u, 0.7),
				u0, opt_new)
end

opts = Cont.ContinuationPar(dsmin = 0.0001,
			dsmax = 0.0035,
			ds = -0.001,
			doArcLengthScaling = false,
			a = 0.5,
			newtonOptions = opt_new,
			detect_fold = false,
			maxSteps = 2340,
			theta = .5, plot_every_n_steps = 200)
	@assert opts.a<=1.5 "sinon ca peut changer le sens du time step"
	opts.newtonOptions.maxIter = 100
	opts.newtonOptions.tol = 1e-9
	opts.newtonOptions.damped = false
	# opts.detect_fold = true

	# opts.maxSteps = 180
	# opts.theta = .5
	br, u1 = @time Cont.continuation(
					(x, p)->  R_SH(x, p, 2., Lsh),
					(x, p)->Jac_sp(x, p, 2., Lsh),
					sol1, 0.7, opts,
					plot = true,
					plotsolution = (x;kwargs...)->(plot!(X, x, subplot=4, ylabel="solution", label="")))
#####################################################
# case with computation of eigenvalues
# opt_new = Cont.NewtonPar(linsolve = Default(),	eigsolve = eig_KrylovKit{Float64}())
back to top