Untitled
unknown
julia
2 years ago
3.9 kB
4
Indexable
module fDNLS_Direct using DifferentialEquations, GLMakie, LinearAlgebra, OffsetArrays, NonlinearSolve NonlinearSolve.ForwardDiff.can_dual(::Type{ComplexF64}) = true using SpecialFunctions: zeta GLMakie.activate!() #= DifferentialEquations - Used for the ODE Solver GLMakie - 2D/3D plotting backend, LinearAlgebra, SpecialFunctions - Used for importing the Riemann zeta function =# N = 10 # 2N-1 ϵ = 1 p = [ϵ] α = 0.9 ω = 1 γ = 1 tspan = (0.0,60) x=-N+1:N-1 function main(N) #= Per the boundary conditions, solve for u = (u_{-N+1}, ..., u_{N-1}) =# # Per u_(-N) = 0 and u_(N) = 0, then we solve the (2N+1) - 2 = 2N-1 total equations RHS(u,p,t) = -im*(p[1]*coefficientMatrix(α,x)*u - γ*abs2.(u).*u) initial_values = Complex.(ground_state(ω,α,x)).*exp.(im*x) # Equivalent to Ode45 from Matlab prob = ODEProblem(RHS, initial_values, tspan, p) solution_direct = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # At every point in time, ell2 norm of solution. #print(size(solution_direct.u)) #= for i in 1:size(solution_direct.u)[1] println(norm(solution_direct.u[i]), 2) end =# plot_main(solution_direct, x) # Surface and contour plot of |u|² return solution_direct end function coefficientMatrix(α,x) M = OffsetArray(zeros(length(x), length(x)), x,x) # Offset indices for M[i,j]. For example, M[-N+1, -N+1] corresponds to M[1,1]. for i in firstindex(M,1) : lastindex(M,1) for j in axes(M,2) i == j ? M[i,i] = 1/(abs(i - (-N))^(1+α)) + sum(1/(abs(i - m)^(1+α)) for m in x if m!=i) + 1/(abs(i - N)^(1+α)) : M[i,j] = -1/(abs(i - j)^(1+α)) end end return M.parent end function onSite(ω,α,x) # Assumption: qₙ = q₋ₙ, q₀ >> 1 >> q₁ >> ... Q = OffsetVector(zeros(length(x)),x) Q[0] = sqrt(ω + 2*zeta(1+α)) Q[1] = Q[0]/(2*zeta(1+α) - (2)^(-1-α) + ω) Q[-1] = Q[1] # Obtain q₂, q₃, ... via equation (3.1) for n in 2:maximum(x) #drop first 1 since already initialized #asymptotic_onsite = sum(Q[j]/((n+j)^(1+α)) for j in 1-n:n-1)\ asymptotic_onsite = (Q[0]/(n^(1+α)) + sum((1/abs(n-j)^(1+α) + 1/abs(n+j)^(1+α))*Q[j] for j in 1:n-1))/(2*zeta(1+α) - (2n)^(-1-α) + ω) Q[-n] = asymptotic_onsite Q[n] = asymptotic_onsite end return Float64.(Q.parent) # typeof(Q) = OffsetVector, typeof(Q.parent) = Vector end function ground_state(ω, α, x) # 0 = ω*qₙ + (Lq)ₙ - qₙ³ p = [ω] RHS(u,p) = p[1]*u + coefficientMatrix(α,x)*u - u.^3 initialValues = onSite(p[1],α,x) # retrieves the asymptoic onsite sequence Q = [q₋ₙ, ..., q₋₂, q₋₁, q₀, q₁, ..., qₙ] probN = NonlinearProblem(RHS, initialValues, p) sol = solve(probN, NewtonRaphson(), reltol = 1e-9) return sol end function plot_main(solution_direct,x) fig = Figure(backgroundcolor=:snow2) y = solution_direct.t z= [abs(solution_direct[i,j])^2 for i in axes(solution_direct,1), j in axes(solution_direct,2)] # solution_direct[i,j] is the ith component at timestep j. surface_axis=Axis3(fig[1:2,1],xlabel="Space", ylabel="Time", zlabel="|u|²", title="Intensity") p1 = surface!(surface_axis, x, y, z) axis_initial_condition = Axis(fig[1,2],xlabel="space",title="Initial condition: ''fsolve'' with onsite initial", subtitle=L"-\omega q_{n} = \epsilon\sum_{m \in \mathbb{Z}, m\neq n}{\frac{q_n - q_m}{|n-m|^{1+\alpha}}} - |q_n|^2 q_n") p2 = scatter!(axis_initial_condition, x, ground_state(ω,α,x)) norm_of_solution = [norm(solution_direct.u[i],2) for i in 1:size(solution_direct.u)[1]] axis_of_norm = Axis(fig[2,2],xlabel="Time",title="Conservation of 2-norm of solution") p3 = scatter!(axis_of_norm, solution_direct.t, norm_of_solution) display(fig) end #const N=5 main(N) end
Editor is loading...