How to avoid issue of unstable solution in case of Differential Equation based models

New to Julia
,
1

Hello,

I was curious that how to approach this problem with unstable solution. I have already tried different solver such as TRBDF2(), Rodas5() , Rosenbrock23() but it didn’t worked.

using DifferentialEquations 
using NaNMath; nm=NaNMath


function ab_model(du, u, p, t) 

    k_1, k_2, n     = p
    A, B, ab        = u     
  
    nA = 1e2

    du[1]   = -(k_1 * A * nA * max(0,nm.pow((B / nA), n)) )+ (k_2 * ab)          
    du[2]   = -(k_1 * A * nA * max(0,nm.pow((B / nA), n))) + (k_2 * ab)         
    du[3]   =  (k_1 * A * nA * max(0,nm.pow((B / nA), n))) - (k_2 * ab)  

  
end



A_ = [1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3]
u0 = [A_[6], 1e-2, 1e-6]
p =  [1, 0.00001, 0.07]

prob  = ODEProblem(ab_model ,u0 , (0.0, 10), p )
sol  = solve(prob, Rodas4(), saveat= 1)

┌ Warning: At t=2.0489651529636236e-6, dt was forced below floating point epsilon 4.235164736271502e-22, and step error estimate = NaN. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase ~/.julia/packages/SciMLBase/7BHQj/src/integrator_interface.jl:623
retcode: Unstable
Interpolation: 1st order linear
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [100.0, 0.01, 1.0e-6]
2

have you plotted to see what the solution looks like? the issue may be mathematical. plenty of ODEs blow up to infinity in finite time.

3

my test problem mentioned in above examples is solved with use of BigFloat or keeping power expression out of equation. Though in actual model it still doesn’t help.

# Define our model function with improved stability
function ab_model(du, u, p, t) 
    k_1, k_2, n = p  # Unpack parameters
    A, B, ab = u     # Unpack state variables
  
    nA = 1e2  # Constant
    
    # Calculate the power term safely
    # Using max to prevent negative numbers being raised to non-integer powers
    power_term = max(0, B / nA)^n
    
    # Define the differential equations
    du[1] = -k_1 * A * nA * power_term + (k_2 * ab)          
    du[2] = -k_1 * A * nA * power_term + (k_2 * ab)          
    du[3] = k_1 * A * nA * power_term - (k_2 * ab)  
end
4

Can you show the bigfloat case that works?

5

I just checked this is still causing unstability. except for simple initial condition such as u = [1, 1, 0] while unstable for u = [1e2, 1e-2, 1e5] with same parameters p = [1, 1e-5, 0.07]

function ab_model(du, u, p, t) 
    k_1, k_2, n = p  # Unpack parameters
    A, B, ab = u     # Unpack state variables
  
    nA = BigFloat(1e2)  # Using BigFloat for all calculations
    
    # Convert all values to BigFloat for consistent precision
    A_big = BigFloat(A)
    B_big = BigFloat(B)
    ab_big = BigFloat(ab)
    k1_big = BigFloat(k_1)
    k2_big = BigFloat(k_2)
    n_big = BigFloat(n)
    
    # # Calculate power term safely
    # power_term = (B_big / nA)^n_big
    # if isnan(power_term) || power_term < 0
    #     power_term = BigFloat(0.0)
    # end
    
    # Define the differential equations with BigFloat calculations
    du[1] = -k1_big * A_big * nA * (B_big / nA)^n_big + (k2_big * ab_big)          
    du[2] = -k1_big * A_big * nA * (B_big / nA)^n_big + (k2_big * ab_big)          
    du[3] = k1_big * A_big * nA * (B_big / nA)^n_big - (k2_big * ab_big)
    
    # Convert back to regular floats for the ODE solver
    du[1] = Float64(du[1])
    du[2] = Float64(du[2])
    du[3] = Float64(du[3])
end
6

Oh so the solver isn’t set to use BigFloats, just an intermediate computation? Yeah… that’s just a big ^ :sweat_smile: