# Why won't this simple differential equation work correctly?

**URL:** <https://discourse.julialang.org/t/why-wont-this-simple-differential-equation-work-correctly/132344>\
**Category:** General Usage\
**Tags:** question, differentialequation\
**Created:** [September 13, 2025, 11:59am UTC](https://discourse.julialang.org/t/why-wont-this-simple-differential-equation-work-correctly/132344 "2025-09-13T11:59:39Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [September 13, 2025, 11:59am UTC](https://discourse.julialang.org/t/why-wont-this-simple-differential-equation-work-correctly/132344/1 "2025-09-13T11:59:39Z")

</div>

Hello,

I am new to differential equations and I am trying to implement as simple differential equation model describe [here](https://academic.oup.com/pnasnexus/article/4/3/pgaf082/8069205). 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**
>
> ```julia-auto
> 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**
>
> ```julia-auto
> 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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 13, 2025, 12:23pm UTC](https://discourse.julialang.org/t/why-wont-this-simple-differential-equation-work-correctly/132344/2 "2025-09-13T12:23:02Z")

</div>

> [@Christopher\_Fisher](#):
>
> I have no idea what this normal cdf is. Why is it there?

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](https://en.wikipedia.org/wiki/Sigmoid_function).

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](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) 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.)

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [September 13, 2025, 1:03pm UTC](https://discourse.julialang.org/t/why-wont-this-simple-differential-equation-work-correctly/132344/3 "2025-09-13T13:03:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 13, 2025, 3:27pm UTC](https://discourse.julialang.org/t/why-wont-this-simple-differential-equation-work-correctly/132344/4 "2025-09-13T15:27:18Z")

</div>

> [@Christopher\_Fisher](#):
>
> because DiffEQ handles the continuous approximation?

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.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [September 13, 2025, 4:42pm UTC](https://discourse.julialang.org/t/why-wont-this-simple-differential-equation-work-correctly/132344/5 "2025-09-13T16:42:11Z")

</div>

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:

```julia-auto
    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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 13, 2025, 7:30pm UTC](https://discourse.julialang.org/t/why-wont-this-simple-differential-equation-work-correctly/132344/6 "2025-09-13T19:30:52Z")

</div>

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.