# How to avoid issue of unstable solution in case of Differential Equation based models
**URL:** https://discourse.julialang.org/t/how-to-avoid-issue-of-unstable-solution-in-case-of-differential-equation-based-models/128714
**Category:** New to Julia
**Tags:** question, differentialequation
**Created:** [May 5, 2025, 10:14am UTC](https://discourse.julialang.org/t/how-to-avoid-issue-of-unstable-solution-in-case-of-differential-equation-based-models/128714 "2025-05-05T10:14:59Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@Sahil\_Khan](https://discourse.julialang.org/u/Sahil_Khan)
#### Post date: [May 5, 2025, 10:14am UTC](https://discourse.julialang.org/t/how-to-avoid-issue-of-unstable-solution-in-case-of-differential-equation-based-models/128714/1 "2025-05-05T10:14:59Z")
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.
```julia
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)
```
```julia
┌ 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]
```
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [May 5, 2025, 11:57am UTC](https://discourse.julialang.org/t/how-to-avoid-issue-of-unstable-solution-in-case-of-differential-equation-based-models/128714/2 "2025-05-05T11:57:37Z")
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.
---
### Author:  [@Sahil\_Khan](https://discourse.julialang.org/u/Sahil_Khan)
#### Post date: [May 5, 2025, 1:20pm UTC](https://discourse.julialang.org/t/how-to-avoid-issue-of-unstable-solution-in-case-of-differential-equation-based-models/128714/3 "2025-05-05T13:20:51Z")
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.
```julia
# 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
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 5, 2025, 1:27pm UTC](https://discourse.julialang.org/t/how-to-avoid-issue-of-unstable-solution-in-case-of-differential-equation-based-models/128714/4 "2025-05-05T13:27:48Z")
> [@Sahil\_Khan](#):
>
> t 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 h
Can you show the bigfloat case that works?
---
### Author:  [@Sahil\_Khan](https://discourse.julialang.org/u/Sahil_Khan)
#### Post date: [May 5, 2025, 1:38pm UTC](https://discourse.julialang.org/t/how-to-avoid-issue-of-unstable-solution-in-case-of-differential-equation-based-models/128714/5 "2025-05-05T13:38:46Z")
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]`
```julia
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
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 5, 2025, 4:33pm UTC](https://discourse.julialang.org/t/how-to-avoid-issue-of-unstable-solution-in-case-of-differential-equation-based-models/128714/6 "2025-05-05T16:33:22Z")
Oh so the solver isn’t set to use `BigFloat`s, just an intermediate computation? Yeah… that’s just a big `^` 😅