Simulating a system of ODES in Julia with ODEProblem

New to Julia
, ,
1

I am using ODEProblem to solve a system of ODEs: https://mathb.in/75542 such that u_1 = u_N = 0.

To do this, I wrote the code shown at the end of the post to generate this solution: https://i.imgur.com/wQC8aTj.png I expected the solution to decay near the boundary, but instead we have quite the opposite.

**Question: ** Does the following code indeed solve the system (1) shown in https://mathb.in/75542 ?


using DifferentialEquations, GLMakie, LinearAlgebra, SpecialFunctions
GLMakie.activate!()

#=
DifferentialEquations - Used for the ODE Solver
GLMakie               - 2D/3D plotting backend,
LinearAlgebra,
SpecialFunctions      - Used for importing the Riemann zeta function
=#

const N = 31    

centerGrid  = (N-1)/2;
x = -centerGrid:centerGrid;

#f(n,α) = sum(g(n,m,α) for m in 1:N if m!=n) # Diagonal elements for the coefficient matrix M
#g(n,m,α) = 1/abs(n-m)^(1+α) # Off-diagonal elements for the coefficient matrix M

function main()

    p   = [1] # ODEProblem() requries a parameter list: `p`. Here, p = [h]
    α   = 0.2
    h   = 1
   
    tspan = (0.0,100)

    RHS(u,p,t) = 1*im*(p[1]*coefficientMatrix(α)*u - @.(((abs(u))^2).*u) )

    initial_values = Complex.(sech.(x))

    # Equivalent to Ode45 from Matlab
    prob = ODEProblem(RHS, initial_values, tspan, [h])
    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)   

    plot_main(sol) # Surface and contour plot of |u|²
    return sol
end

function diagonal_coefficients(n,α)
    if n == 1 || n == N
        diag_value = 0
    else
        diag_value = 1/abs(n-1)^(1+α) + sum(1/abs(n-m)^(1+α) for m in 2:N-1 if m!=n) + 1/abs(n-N)^(1+α) # Entry for M[n,n]. row n, col n
    end
        return diag_value
end

function off_diagonal_coefficients(n,m,α)
    if n == 1 || n == N
        off_diag_value = 0
    else
        off_diag_value = -1/abs(n-m)^(1+α) # Entry for M[n,m]. row n, col m
    end
    return off_diag_value
end

function coefficientMatrix(α) 
    # M is the coefficient matrix corresponding to the 1-dimensional nonlinear discrete operator of the fDNLS
    M = Matrix{Float64}(undef, length(x), length(x))
    for i in 1:size(M,1) # Iterate through the rows
        for j in 1:size(M,2) # Iterate through the columns
            if i==j
                M[i,i] = diagonal_coefficients(i, α)
            elseif i!=j
                M[i,j] = off_diagonal_coefficients(i,j,α)
            end
        end
    end

    return M
end

function plot_main(sol)
    fig = Figure(backgroundcolor=:snow2)

    y = sol.t
    z= [abs(sol[i,j])^2  for i=1:size(sol)[1], j=1:size(sol)[2]]  # sol[i,j] is the ith component at timestep j. # sol[i,j] is the ith component at timestep j. 

    p1 =  surface(fig[1,1], x, y, z,axis=(type=Axis3, xlabel="Space", ylabel="Time", zlabel="|u|²", xticks=-centerGrid:2:centerGrid))
    p2 =  contourf(fig[1,2], x, y, z, axis=(xlabel="Space", ylabel="Time"))
    display(fig)
end

main()
2

@xtalax can you take a stab here? If not I set a reminder for Thursday

3

Update: Well, the code works but now I don’t understand the math behind why we are only solving for N-2 equations. The problem with the above code were the initial conditions used:

initial_values = Complex.(sech.(x))

does not work well and causes instability, but different initial conditions (which work somewhat by construction and whose explaination are outside the scope of this post) work fine:

     initial_values = Complex.(onSite(30,α)) # onSite(ω, α)

Below you can examine the updated code:

using DifferentialEquations, GLMakie, LinearAlgebra, SpecialFunctions
GLMakie.activate!()

#=
The fDNLS is written as iuₙ' = ∑(uₙ - uₘ)/|n-m|^(1+α) - |uₙ|^2 uₙ, where the sum is evaluated such that m \in 1:N and m != n. 
Dirichlet boundary values may be imposed: u₁ = 0, uₙ = 0.

Standing wave solutions are of the form u(x,t) = q(x)exp[i ω t], where q(x) here is the "Ground State Solution" obtained by solution by the 
corresponding system of nonlinear equations after substituting the ansatz. To solve the above system of ODEs, take u(x,0) = q(x). Ref. https://royalsocietypublishing.org/doi/epdf/10.1098/rspa.2014.0364
=#

#=
DifferentialEquations - Used for the ODE Solver
GLMakie               - 2D/3D plotting backend,
LinearAlgebra,
SpecialFunctions      - Used for importing the Riemann zeta function
=#

const N = 31    

centerGrid  = (N-1)/2;
x = -centerGrid:centerGrid;

h(j,n,α) = 1/(n-j)^(1+α) + 1/(n+j)^(1+α) # Used in the onSite() function
Q₊  = zeros(Int(centerGrid+1))

function main()

    p   = [1] # ODEProblem() requries a parameter list: `p`. Here, p = [h]
    α   = 0.2
    h   = 1
   
    tspan = (0.0,10.0)
    A = Complex.(diagm(ones(length(x))))
    A[1,1] = 0
    A[length(x),length(x)] = 0 

    RHS(u,p,t) = 1*im*(p[1]*coefficientMatrix(α)*u - A*@.(((abs(u))^2)*u) )

    #initial_values = Complex.(sech.(x))
    initial_values = Complex.(onSite(30,α)) # onSite(ω, α)
    # Equivalent to Ode45 from Matlab
    prob = ODEProblem(RHS, initial_values, tspan, [h])
    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)   

    plot_main(sol) # Surface and contour plot of |u|²
    return sol
end

function diagonal_coefficients(n,α)
    
    if n == 1 || n == N
        diag_value = 0
    else
        diag_value = 1/abs(n-1)^(1+α) + sum(1/abs(n-m)^(1+α) for m in 2:N-1 if m!=n) + 1/abs(n-N)^(1+α) # Entry for M[n,n]. row n, col n
    end
    
    return diag_value
end

function off_diagonal_coefficients(n,m,α)
    
    if n == 1 || n == N
        off_diag_value = 0
    else
        off_diag_value = 1/abs(n-m)^(1+α) # Entry for M[n,m]. row n, col m
    end
    
    return off_diag_value
end

function coefficientMatrix(α) 
    # M is the coefficient matrix corresponding to the 1-dimensional nonlinear discrete operator of the fDNLS
    M = Matrix{Float64}(undef, length(x), length(x))
    for i in 1:size(M,1) # Iterate through the rows
        for j in 1:size(M,2) # Iterate through the columns
            if i==j
                M[i,i] = diagonal_coefficients(i, α)
            elseif i!=j
                M[i,j] = off_diagonal_coefficients(i,j,α)
            end
        end
    end

    return M
end

function onSite(ωᵢ,α)
    # Assumption: qₙ = q₋ₙ, q₀ >> 1 >> q₁ >> ...
    q₀ = sqrt(ωᵢ + 2*zeta(1+α))
    q₁ = q₀/(2*zeta(1+α) - 1/2^(1+α) + ωᵢ)

    Q₊[1] = q₀
    Q₊[2] = q₁

    # Obtain q₂, q₃, ... via equation (3.1)
    for n in 3:size(Q₊)[1]
        Q₊[n] = (q₀/(n)^(1+α) + sum(h(j,n,α)*Q₊[j+1] for j in 1:n-1))/(2*zeta(1+α) - 1/(2*n)^(1+α) + ωᵢ)
    end

    Q₋ = reverse(Q₊[2:size(Q₊)[1]]) # [q₋ₙ, ..., q₂, q₁]. Does not include q₀. 
    Q = append!(Q₋, Q₊)             # [q₋ₙ, ..., q₂, q₁, q₀, q₁, ..., qₙ] 

    return Q
end

function plot_main(sol)
    fig = Figure(backgroundcolor=:snow2)

    y = sol.t
    z= [abs(sol[i,j])^2  for i=1:size(sol)[1], j=1:size(sol)[2]]  # sol[i,j] is the ith component at timestep j. # sol[i,j] is the ith component at timestep j. 

    p1 = surface(fig[1,1], x, y, z,axis=(type=Axis3, xlabel="Space", ylabel="Time", zlabel="|u|²", xticks=-centerGrid:2:centerGrid, title="Intensity plot with u(x,0) = q(x)"))
    p2 =  contourf(fig[1,2], x, y, z, axis=(xlabel="Space", ylabel="Time"))
    display(fig)
end

main()

My math question is here: Why do Dirichlet boundary conditions imply we solve these $N-2$ equations instead of $N$ equations? - Mathematics Stack Exchange

4

Mathematically the PDE is defined at the interior and the boundary conditions give explicit relationships for the 2 boundary variables, so you only have N-2 equations, one for each of the non-boundary terms.

5

Thanks. I guess what I do in the code works, but its pretty messy. I’d appreciate any cleaning up if see something