Why won't this simple differential equation work correctly?

General Usage
,
1

Hello,

I am new to differential equations and I am trying to implement as simple differential equation model describe here. The basic model is described in equation 4. I am able to reproduce the expected behavior of case 3 in Figure 1 when using the parameters. r = .20, δ = 0.00, α = .10, β = 0.20. However, I am not able to match the behavior when using different parameters that satisfy case 3, as shown in the Julia code. I expected the curves to have the same qualitative properties.

Julia Code 
using ComponentArrays
import DifferentialEquations as DE
using Plots
using Test

function ode_model!(
    du::AbstractVector{<:R},
    u::AbstractVector{<:R},
    θ::ComponentArray{<:R},
    _
) where {R <: Real}
    (; δ, r, α, β) = θ
    p1 = (α * (1 - r) * (2 * u[1] - 1) - β * r * (2 * u[2] - 1)) > δ
    p2 = (α * (1 - r) * (2 * u[1] - 1) - β * r * (2 * u[2] - 1)) < -δ
    p3 = (α * (2 * u[2] - 1) - β * (1 - r) * (2 * u[1] - 1)) > δ
    p4 = (α * (2 * u[2] - 1) - β * (1 - r) * (2 * u[1] - 1)) < -δ

    # change in behavior for blue node
    du[1] = (1 - u[1]) * p1 - u[1] * p2
    # change in behavior for red node
    du[2] = (1 - u[2]) * p3 - u[2] * p4
    return nothing
end

r = .37
δ = 0.00
α = .74
β = 0.92
# ensure parameters satisfy conditions for case 3
@test (r / (1 - r)) < (β / α)
@test (r / (1 - r)) < (α / β)
θ = ComponentArray(; r, δ, α, β)
u0 = [0.80, 0.80]
tspan = (0.0, 10.0)
prob = DE.ODEProblem(ode_model!, u0, tspan, θ)
sol = DE.solve(prob)
plot(sol; idxs = 1, ylims = (0, 1))
plot!(sol; idxs = 2, linestyle = :dash)

My Julia example does not match the results of the code Python associated with the paper. I extracted the necessary components below, and disabled the parameter rho because its not relevant for the basic model.

Python Code 
def Simulate_AP_model_SBM_Graph(Initial_H1_prob_blue, Initial_H1_prob_red, Alpha_b, Beta_b, Rho_b, Intertia_b, Alpha_r, Beta_r, Rho_r, Intertia_r, R, n = 5000000, T = 10, delta = 0.001):
    
    theta_t_blue = [Initial_H1_prob_blue]
    theta_t_red = [Initial_H1_prob_red]

    T = 10
    time = np.arange(0,T,delta)
    for t in time[:-1]:
        (theta_t_blue_dot, theta_t_red_dot) = AP_model_SBM_Graph(theta_t_blue[-1],theta_t_red[-1], t, Alpha_b, Beta_b, Rho_b, Inertia_b, Alpha_r, Beta_r, Rho_r, Inertia_r, R, n)
        theta_t_blue.append(theta_t_blue[-1] + delta*theta_t_blue_dot)
        theta_t_red.append(theta_t_red[-1] + delta*theta_t_red_dot)
    
    return time, theta_t_blue, theta_t_red

def AP_model_SBM_Graph(theta_t_b,theta_t_r, t, Alpha_b, Beta_b, Rho_b, Inertia_b, Alpha_r, Beta_r, Rho_r, Inertia_r, R, n):
    b_in_link_prob = (1-R)#*Rho_b
    b_out_link_prob = R#*(1-Rho_b)
    
    mu_b = np.array([b_in_link_prob*theta_t_b, b_in_link_prob*(1-theta_t_b), b_out_link_prob*theta_t_r, b_out_link_prob*(1-theta_t_r)]).reshape(4,1)
    Sigma_b = np.zeros((4,4))
    for i in [0,1,2,3]:
        Sigma_b[i,i] = mu_b[i,0]*(1-mu_b[i,0])/n

    r_in_link_prob = R#*Rho_r
    r_out_link_prob = (1-R)#*(1-Rho_r)                         

    mu_r = np.array([r_in_link_prob*theta_t_r, r_in_link_prob*(1-theta_t_r), r_out_link_prob*theta_t_b, r_out_link_prob*(1-theta_t_b)]).reshape(4,1)
    Sigma_r = np.zeros((4,4))
    for i in [0,1,2,3]:
        Sigma_r[i,i] = mu_r[i,0]*(1-mu_r[i,0])/n
                                        
                    
    coeffs_blue = np.array([Alpha_b, -Alpha_b, -Beta_b, Beta_b]).reshape(4,1)
    coeffs_red = np.array([Alpha_r, -Alpha_r, -Beta_r, Beta_r]).reshape(4,1)    
                    
    p_b_01 = 1-NormalDist(mu=np.matmul(coeffs_blue.transpose(),mu_b)-Inertia_b, sigma=np.sqrt(np.matmul(np.matmul(coeffs_blue.transpose(),Sigma_b),coeffs_blue))).cdf(0)
    p_b_10 = NormalDist(mu=np.matmul(coeffs_blue.transpose(),mu_b)+Inertia_b, sigma=np.sqrt(np.matmul(np.matmul(coeffs_blue.transpose(),Sigma_b),coeffs_blue))).cdf(0)                          
        
    p_r_01 = 1-NormalDist(mu=np.matmul(coeffs_red.transpose(),mu_r)-Inertia_r, sigma=np.sqrt(np.matmul(np.matmul(coeffs_red.transpose(),Sigma_r),coeffs_red))).cdf(0)
    p_r_10 = NormalDist(mu=np.matmul(coeffs_red.transpose(),mu_r)+Inertia_r, sigma=np.sqrt(np.matmul(np.matmul(coeffs_red.transpose(),Sigma_r),coeffs_red))).cdf(0)  
    theta_t_blue_dot = (1-theta_t_b)*p_b_01 - theta_t_b*p_b_10
    theta_t_red_dot = (1-theta_t_r)*p_r_01 - theta_t_r*p_r_10                        
                                                
    return [theta_t_blue_dot, theta_t_red_dot]   

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.stats import bernoulli
from scipy.stats import pareto
from statistics import NormalDist
import networkx as nx
import matplotlib.pyplot as plt
import random
import math

# Set initial conditions
Initial_H1_prob_blue = 0.8
Initial_H1_prob_red = 0.8
# Set parameters
R = 0.37
Alpha = 0.74
Beta = 0.92

Rho = .5
Inertia = 0.0

Alpha_b = Alpha
Beta_b = Beta
Rho_b = Rho
Inertia_b = Inertia
Alpha_r = Alpha
Beta_r = Beta
Rho_r = Rho
Inertia_r = Inertia

time, theta_t_blue, theta_t_red = Simulate_AP_model_SBM_Graph(Initial_H1_prob_blue, Initial_H1_prob_red, Alpha_b, Beta_b, Rho_b, Inertia_b, Alpha_r, Beta_r, Rho_r, Inertia_r, R, n = 5000000, T = 10, delta = 0.001)
plt.plot(time, theta_t_blue)
plt.plot(time, theta_t_red)

I noticed two main differences. For some reason the equations include a normal cdf, which is not described in paper (see for example p_b_01 in AP_model_SBM_Graph), and the model in Python is less sensitive to changes in parameter values (perhaps due to the use of the normal cdf).

I have no idea what this normal cdf is. Why is it there? Did I set up the code incorrectly, or am I using the wrong type of solver? Any help would be appreciated.

2

Your equations have step-function discontinuities, which are difficult for standard ODE solvers to handle (the algorithms assume smooth functions). It looks like the normal CDF replaces this discontinuity with a smooth sigmoid approximation.

You can, of course, do the same thing in Julia — replace NormalDist(mu=x, sigma=y).cdf(0) with StatsFuns.normcdf(x, y, 0) from the StatsFuns.jl package.

If you want to handle the exact discontinuities in Julia, you will typically want to do so by one or more continuous callbacks that tell the solver where the discontinuities occur so that it can handle them efficiently (and with high-order accuracy).

(It doesn’t help that the Python code is written in an extremely convoluted way, in a desperate and probably misguided attempt to “vectorize” all of the operations, even though the resulting matrices are so tiny that the overhead of creating them probably swamps the gains of using numpy. So it will require some patience to analyze exactly what it is doing.)

4 Likes
3

Thank you for your reply. I had a suspicion that the step functions were part of the problem, but I wasn’t quite sure.

I agree that the Python code is difficult to parse. I’ll see if I can get one or both of the methods to work. The main thing I don’t understand about the sigmoid approximation is how the sigma parameter should be calculated in DifferentialEquations.jl. Do you happen to know if I can drop the n in Sigma_r[i,i] = mu_r[i,0]*(1-mu_r[i,0])/n because DiffEQ handles the continuous approximation? It seems like selecting the right value is important and I wonder if it is too small in the Python code, because the model predictions seem somewhat insensitive to changes in the parameters.

4

DiffEQ doesn’t do the sigmoid approximation for you — it tries to accurately solve the problem you give it. The sigmoid approximation is changing the problem, and that’s up to you.

As I said, DiffEQ should be perfectly capable of solving the “exact” problem with a discontinuous right-hand side, as long as you use continuous callbacks to tell it where the discontinuities occur. (Otherwise, it will struggle to get an accurate result, probably putting a huge number of points at each discontinuity to try to resolve them adaptively.) But if you want to replicate the results of the paper you’ll need to use the same sigmoid relaxation.

1 Like
5

After reading through the documentation, I am interested in using callbacks. It seems like it might be less error prone (or at least easier to setup) than the sigmoid approximation. I think I understand the examples, but I am uncertain how to use it for the current model. If I understand correctly, the problem occurs when either of the following switch from true to false or vice versa:

    p1 = (α * (1 - r) * (2 * u[1] - 1) - β * r * (2 * u[2] - 1)) > δ
    p2 = (α * (1 - r) * (2 * u[1] - 1) - β * r * (2 * u[2] - 1)) < -δ
    p3 = (α * (2 * u[2] - 1) - β * (1 - r) * (2 * u[1] - 1)) > δ
    p4 = (α * (2 * u[2] - 1) - β * (1 - r) * (2 * u[1] - 1)) < -δ

As far as I can tell, there is not a way to evaluate the equations above from condition(u, t, integrator) . integrator has lots of fields, but I am guessing most of the should not be accessed by the user (except u it appears). The other thing that is unclear to me is what affect! should be if anything. Can someone provide any guidance or point me to an example with similar characteristics?

6

You can install two callbacks, one with the condition function (α * (1 - r) * (2 * u[1] - 1) - β * r * (2 * u[2] - 1)) - δ, and the other with the condition (α * (2 * u[2] - 1) - β * (1 - r) * (2 * u[1] - 1)) - δ. These hit 0 when the two discontinuities occur, respectively.

Your affect! function can, for example, swap the values of p1 and p2 (for the first condition) or p3 and p4 (for the second condition), which you can store in integrator.p (the parameters, which you could make a mutable struct, or store [p1, p2, p3, p4] as an array field in a named-tuple parameter or other immutable struct), rather than computing them at each timestep.

1 Like