#NetworkDynamics.jl Example

This example shows the basic usage of AmbientForcing.jl and NetworkDynamics.jl. It is adopted from the Network Dynamics Docs. This example only works with NetworkDynamics v0.5.0 or newer.

We begin with the NetworkDynamics Setup. First we include all necessary packages.

using AmbientForcing
using NetworkDynamics, Distributions
using Graphs
using OrdinaryDiffEq

We generate a random graph with 10 nodes and an average degree of 4.

const N_plastic = 10 # number of nodes
k = 4  # average degree
g = barabasi_albert(N_plastic, k)

{10, 24} undirected simple Int64 graph

The coupling function is modeled by a differential algebraic equation with mass matrix:

\[0 \cdot \frac{de_1}{dt} = e_2 \cdot \frac{sin(v_{s1} - v_{d1} + α)}{N} - e_1\]

which is equivalent to: $e_1 = e_2 \cdot \frac{sin(v_{s1} - v_{d1} + α)}{N}$

This model comes from the following publication: R. Berner et. al., "Multiclusters in Networks of Adaptively Coupled Phase Oscillators." SIAM Journal on Applied Dynamical Systems 18.4 (2019): 2227-2266.

In NetworkDynamics.jl we define the edge dynamics as following:

function kuramoto_plastic_edge!(de, e, v_s, v_d, p, t)
    de[1] =  e[2] * sin(v_s[1] - v_d[1] + α) / N_plastic - e[1]
    de[2] = - ϵ * (sin(v_s[1] - v_d[1] + β) + e[2])

    nothing
end

kuramoto_plastic_edge! (generic function with 1 method)

and the node dynamics as:

function kuramoto_plastic_vertex!(dv, v, edges, p, t)
    dv .= 0
    for e in edges
        dv .-= e[1]
    end
end

kuramoto_plastic_vertex! (generic function with 1 method)

We define the following constants:

const ϵ = 0.1
const α = .2π
const β = -.95π

-2.9845130209103035

Then we generate a network dynamics problem:

plasticvertex = ODEVertex(f = kuramoto_plastic_vertex!, dim =1)
mass_matrix_plasticedge = zeros(2,2)
mass_matrix_plasticedge[2,2] = 1. # First variables is set to 0

plasticedge = ODEEdge(f = kuramoto_plastic_edge!, dim=2, sym=[:e, :de], coupling=:undirected,mass_matrix = mass_matrix_plasticedge);
kuramoto_plastic! = network_dynamics(plasticvertex, plasticedge, g)

(::SciMLBase.ODEFunction{true, NetworkDynamics.NetworkDE{Graphs.SimpleGraphs.SimpleGraph{Int64}, NetworkDynamics.GraphDataBuffer{SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}, Float64, Float64, Vector{NetworkDynamics.ODEVertex{typeof(Main.kuramoto_plastic_vertex!)}}, Vector{NetworkDynamics.ODEEdge{NetworkDynamics.var"#48#55"{typeof(Main.kuramoto_plastic_edge!), Int64}}}, Nothing}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}) (generic function with 7 methods)

The Ambient Forcing part starts here!

Using a random initial condition x0 violates the constraints! The constraints are fulfilled when $g(x) ≈ 0$.

x0_plastic = rand(106)
g_nd = constraint_equations(kuramoto_plastic!)
sum(g_nd(x0_plastic))

-21.780632108373332

We use zeros as the initial conditions for the ambient forcing algo and access the constraint equation $g$:

x0_plastic = zeros(106)
g_nd = constraint_equations(kuramoto_plastic!)
g_nd(x0_plastic)

48-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

This point fulfills the constraints and we can use it as the initial condition for Ambient Forcing.

We start by perturbing all variables at once:

Frand = random_force(kuramoto_plastic!, [0.0, 1.0], Uniform)
afoprob = ambient_forcing_problem(kuramoto_plastic!, x0_plastic, 2.0, Frand)
z_new = ambient_forcing(afoprob, x0_plastic, 2.0, Frand)

106-element Vector{Float64}:
  0.4881305451801231
  0.661878155821338
  0.9283588557511411
  1.3447213860010103
  1.2871617055039515
  1.3248261590887325
  1.3530887935563112
  0.8197622906010652
  0.7039740744652767
  0.23317744492113407
  ⋮
  1.9275373804955085
  0.19272122309997594
  2.013223608486865
 -0.0033425446543881836
  1.6074037506118066
  0.0173018654314102
  0.17577838211614233
 -0.02732094758908837
  0.5787960170791714

As we can see the constraints are not violated!

sum(g_nd(z_new))

8.468113254947278e-8

Now we only perturb the variables $e_{22}$ and $de_{22}$

idx = idx_exclusive(kuramoto_plastic!, ["e_22", "de_22"])
Frand = random_force(kuramoto_plastic!, [0.0, 1.0], Uniform, idx)
z_new = ambient_forcing(afoprob, x0_plastic, 2.0, Frand)

106-element Vector{Float64}:
  0.0
  0.0
  0.0
  0.0
  0.0
  0.12271057850675962
  0.0
  0.0
  0.0
 -0.12271057850675962
  ⋮
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0

We can see that the constraints are fulfilled!

sum(g_nd(z_new))

1.429707288669313e-8