# How can I solve this 4-equation NonlinearSystem? It gives \`retcode: Unstable\`
**URL:** https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384
**Category:** Modelling & Simulations
**Tags:** modelingtoolkit, nonlinearsolve
**Created:** [September 13, 2024, 8:05pm UTC](https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384 "2024-09-13T20:05:00Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@JohsEnevoldsen](https://discourse.julialang.org/u/JohsEnevoldsen)
#### Post date: [September 13, 2024, 8:05pm UTC](https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384/1 "2024-09-13T20:05:00Z")
I am trying to solve a system of equations with NonlinearSolve.jl and ModelingToolkit.jl.
Specifically a simplified version of Stewart’s model for Acid-Base disturbances in blood ([Independence of Variables in Stewart's Model of the Acid-Base Chemistry of the Blood Plasma - ScienceDirect](https://www.sciencedirect.com/science/article/pii/S1474667015378800)).
The following code returns a meaningless result (a negative H+ concentration), and `retcode: Unstable`.
```julia
using NonlinearSolve
using ModelingToolkit
# constants
Kc = 2.46e-11
Ka = 3e-7
@variables HA A H HCO3
@parameters pCO2 Atot SID
# balance equations
balance_eqs = [
Atot ~ HA + A
SID ~ HCO3 + A
]
# mass action equations
action_eqs = [
H ~ Ka * (HA/A)
H ~ Kc * (pCO2/HCO3)
]
stewart_eqs = [balance_eqs; action_eqs]
@mtkbuild ns_stewart = NonlinearSystem(stewart_eqs)
# Initial values (guesses)
guesses_stewart = (
HCO3 => 24.9,
H => 3.95e-8,
A => 12.1,
HA => 5.4
)
ps = (
Atot => 20,
SID => 20,
pCO2 => 5
)
prob = NonlinearProblem(ns_stewart, guesses_stewart, ps)
sol = solve(prob)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 14, 2024, 7:27am UTC](https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384/2 "2024-09-14T07:27:28Z")
```julia
julia> full_equations(ns_stewart)
1-element Vector{Equation}:
0 ~ (-2.46e-11pCO2) / HCO3 + (3.0e-7(Atot + HCO3 - SID)) / (-HCO3 + SID)
julia> unknowns(ns_stewart)
1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
HCO3
```
From this you can see that the only solution is `HCO3 -> infty`
---
### Author:  [@JohsEnevoldsen](https://discourse.julialang.org/u/JohsEnevoldsen)
#### Post date: [September 14, 2024, 10:45am UTC](https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384/3 "2024-09-14T10:45:24Z")
Thanks! I will try to figure out where the mistake is, and hopefully post a working model.
---
### Author:  [@johhell](https://discourse.julialang.org/u/johhell)
#### Post date: [September 14, 2024, 12:10pm UTC](https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384/4 "2024-09-14T12:10:46Z")
MTK works correct.
Please be more carefully with units of the `parameters` and `initial values` you are using. e.g.
- parameter `Kc` is using `mm-mercury` =\> `1mmHg = 133.322Pa`
- pressure `pCO2 = 5.7kPa` - you used `5Pa` as initial value
- concentrations should also be checked
---
### Author:  [@johhell](https://discourse.julialang.org/u/johhell)
#### Post date: [September 14, 2024, 2:09pm UTC](https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384/5 "2024-09-14T14:09:02Z")
a revised version below:
```julia
using NonlinearSolve
using ModelingToolkit
# !!! used UNITS !!!
# eq/liter
# Pa
# constants
mmHG = 133.322 # Pa
Kc = 2.46e-11/mmHG
Ka = 3e-7
@variables HA A H HCO3
@parameters pCO2 Atot SID
eqs = [
H * A ~ Ka * HA
H * HCO3 ~ Kc * pCO2
Atot ~ HA + A
SID ~ HCO3 + A
]
@mtkbuild ns_stewart = NonlinearSystem(eqs)
display(ns_stewart)
println()
# Initial values (guesses)
guesses_stewart = (
HCO3 => 24.9e-3,
H => 3.95e-8,
A => 17.1e-3,
HA => 2e-3
)
ps = (
Atot => 19e-3,
SID => 41.7e-3,
pCO2 => 5.3e3
)
prob = NonlinearProblem(ns_stewart, guesses_stewart, ps)
sol = solve(prob)
@show sol
@show H = -log10(sol[1])
println("HCO3 = $(sol[2]*1000) meq/liter")
```
NEW Results:
```julia
julia> include("stewart1.jl")
Model ns_stewart with 2 equations
Unknowns (2):
H
HCO3
Parameters (3):
pCO2
Atot
SID
sol = [3.9275240698625016e-8, 0.02489948137605699]
H = -(log10(sol[1])) = 7.405881144702646
HCO3 = 24.89948137605699 meq/liter
```
---
### Author:  [@JohsEnevoldsen](https://discourse.julialang.org/u/JohsEnevoldsen)
#### Post date: [September 16, 2024, 7:56am UTC](https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384/6 "2024-09-16T07:56:08Z")
Thanks! Reordering the equations fixes it. Also, good catch with the units.
However, I still do not understand why this specification works:
```julia
eqs1 = [
H * A ~ Ka * HA
H * HCO3 ~ Kc * pCO2
Atot ~ HA + A
SID ~ HCO3 + A
]
@mtkbuild mod1 = NonlinearSystem(eqs1)
full_equations(mod1)
# 2-element Vector{Equation}:
# 0 ~ 1.8451568383312584e-10pCO2 - H*HCO3
# 0 ~ -Atot - HCO3 + SID + 3.3333333333333335e6H*(-HCO3 + SID)
```
while this does not work:
```julia
eqs2 = [
H ~ Ka * HA / A
H ~ Kc * pCO2 / HCO3
Atot ~ HA + A
SID ~ HCO3 + A
]
@mtkbuild mod2 = NonlinearSystem(eqs2)
full_equations(mod2)
# 1-element Vector{Equation}:
# 0 ~ (1.8451568383312584e-10pCO2) / HCO3 + (-3.0e-7(Atot + HCO3 - SID)) / (-HCO3 + SID)
```
To me, they seem equivalent, but is there something special about how MTK handles `/` ?
---
### Author:  [@johhell](https://discourse.julialang.org/u/johhell)
#### Post date: [September 16, 2024, 8:33am UTC](https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384/7 "2024-09-16T08:33:24Z")
My interpretation:
In `equ1` I used the equations I found in the article you referred. The variable `H` here is a part of the LHS expression.
In `equ2` the variable `H` is explicit mentioned in 2 equations. For that reason `H` was eliminated by MTK.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 16, 2024, 11:23am UTC](https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384/8 "2024-09-16T11:23:32Z")
In one case H is nonlinearly related to another variable. It’s just a case the tearing algorithm does not handle yet. In fact, it’s hard to reduce because it’s not generally valid to divide like that because in theory that unknown could be zero.
---
### Author:  [@JohsEnevoldsen](https://discourse.julialang.org/u/JohsEnevoldsen)
#### Post date: [September 17, 2024, 12:29pm UTC](https://discourse.julialang.org/t/how-can-i-solve-this-4-equation-nonlinearsystem-it-gives-retcode-unstable/119384/9 "2024-09-17T12:29:52Z")
Thanks for the explanation. I hope I’ll be able to recognize it if I run into a similar problem in the future.