# Solving System of nonlinear Equations using ModelingToolkit.jl and NonlinearSolve.jl

**URL:** <https://discourse.julialang.org/t/solving-system-of-nonlinear-equations-using-modelingtoolkit-jl-and-nonlinearsolve-jl/59096>\
**Category:** Numerics\
**Created:** [April 12, 2021, 11:43am UTC](https://discourse.julialang.org/t/solving-system-of-nonlinear-equations-using-modelingtoolkit-jl-and-nonlinearsolve-jl/59096 "2021-04-12T11:43:01Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![wst](https://avatars.discourse-cdn.com/v4/letter/w/e79b87/32.png) [@wst](https://discourse.julialang.org/u/wst)\
**Post date:** [April 12, 2021, 11:43am UTC](https://discourse.julialang.org/t/solving-system-of-nonlinear-equations-using-modelingtoolkit-jl-and-nonlinearsolve-jl/59096/1 "2021-04-12T11:43:01Z")

</div>

Hey there,  
I’m fairly new to Julia and this forum (in fact my first post here).  
About my problem:  
I’m currently trying to solve a system of 3 nonlinear equations that emerged as the solution of the hanging chain problem (following [this](https://www.youtube.com/watch?v=EHKFbl3VwMo) video).  
My studies didn’t involve solving NL systems very much, so I naively defined the system following the MTK documentation. But when I try to solve it I’m only getting a vector of 3 NaNs as solution.  
If someone could help and point out my error I would really appreciate it 🙂

In case my versions are relevant:

- Julia: 1.6.0
- ModelingToolkit: 5.14.0
- NonlinearSolve: 0.3.8

Here is my code:

```julia
using ModelingToolkit, NonlinearSolve

@variables x y c_1 c_2 λ
@parameters x_1 x_2 y_1 y_2 m g L_0

f(x, y) = c_1/(m*g) * cosh( m*g/c_1*( x + c_2 ) ) - λ/c_1 - y
n = c_1/(m*g) * ( sinh(m*g/c_1*x_2 + c_2) - sinh(m*g/c_1 * x_1 + c_2) ) - L_0

f_1 = f(x_1, y_1)
f_2 = f(x_2, y_2)

eqs = [
        0 ~ f_1
        0 ~ f_2
        0 ~ n
        ]

ns = NonlinearSystem(eqs, [c_1, c_2, λ], [x_1, x_2, y_1, y_2, m, g, L_0])
guess = [
            c_1 => 1.0
            c_2 => 1.0
            λ => 1.0
        ]

ps = [
        x_1 => 0.0
        x_2 => 500.0e3
        y_1 => 0.0
        y_2 => 0.0
        m => 257.0 * 2.7
        g => 9.81e3
        L_0 => 600.0e3
        ]

prob = NonlinearProblem(ns, guess, ps)
sol = solve(prob, NewtonRaphson())

```

(Edit before posting this:)  
After writing this post it came to my mind, that maybe the System is ill-conditioned with the given parameters. So I set all my current non-zero parameters to 1.0 (L\_0 to 1.1 as its the length of the chain) and indeed I got a non-NaN result. Could someone more proficient in the numerics of nonlinear systems point me in the right direction in setting this up please?

Thanks for your time!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 12, 2021, 2:37pm UTC](https://discourse.julialang.org/t/solving-system-of-nonlinear-equations-using-modelingtoolkit-jl-and-nonlinearsolve-jl/59096/2 "2021-04-12T14:37:20Z")

</div>

This is a @YingboMa question.

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [April 12, 2021, 3:07pm UTC](https://discourse.julialang.org/t/solving-system-of-nonlinear-equations-using-modelingtoolkit-jl-and-nonlinearsolve-jl/59096/3 "2021-04-12T15:07:28Z")

</div>

Did you get the equations right? I recrafted to using NLSolve on a plain Julia function and also got immediate Infs (see below), but then started looking at the equations, and they don’t make sense dimensionally.

E.g. you have `m*g/c_1*( x + c_2 )` as the argument of a cosh and `m*g/c_1*x_2 + c_2` as the argument of a sinh. `c_2` would have to have dimension of `x` (length) in the first but be nondimensional in the second. Are the parentheses misplaced or missing?

```julia
julia> x1, x2, y1, y2 = 0, 500e03, 0, 0
(0, 500000.0, 0, 0)

julia> m, g, L0 = 257*2.7, 9.81e03, 600e03
(693.9000000000001, 9810.0, 600000.0)

julia> f(c) = [c[1]/(m*g) * cosh(m*g/c[1]*(x1 + c[2])) - c[3]/c[1] - y1; 
               c[1]/(m*g) * cosh(m*g/c[1]*(x2 + c[2])) - c[3]/c[1] - y2; 
               c[1]/(m*g) * (sinh(m*g/c[1]*x2 + c[2]) - sinh(m*g/c[2] * x1 + c[2])) - L0]

f (generic function with 1 method)

julia> using NLsolve

julia> cguess = [1.0; 1.0; 1.0]
3-element Vector{Float64}:
 1.0
 1.0
 1.0

julia> solution = nlsolve(f, cguess)
ERROR: During the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number: [1, 2, 3]
Stacktrace:
 [1] check_isfinite(x::Vector{Float64})
   @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/utils.jl:24
 [2] trust_region_(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Vector{Float64}})
   @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:121
 [3] trust_region (repeats 2 times)
   @ ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:235 [inlined]
 [4] nlsolve(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
   @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:26
 [5] nlsolve(f::Function, initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:52
 [6] nlsolve(f::Function, initial_x::Vector{Float64})
   @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:46
 [7] top-level scope
   @ REPL[6]:1

julia> f(cguess)
3-element Vector{Float64}:
 Inf
 Inf
 Inf

julia> cguess = [.1; .1; 0.1]
3-element Vector{Float64}:
 0.1
 0.1
 0.1

julia> f(cguess)
3-element Vector{Float64}:
 Inf
 Inf
 Inf

```

---

<div class="post-metadata">

**Author:** ![wst](https://avatars.discourse-cdn.com/v4/letter/w/e79b87/32.png) [@wst](https://discourse.julialang.org/u/wst)\
**Post date:** [April 12, 2021, 3:26pm UTC](https://discourse.julialang.org/t/solving-system-of-nonlinear-equations-using-modelingtoolkit-jl-and-nonlinearsolve-jl/59096/4 "2021-04-12T15:26:22Z")

</div>

Thanks, you are right the dimensions of c\_2 are wrong.

I think in my function n the term m\*g/c\_1 (=constant) got absorbed by c\_2 and in function f I forgot to do this. I’ll try to fix this and report back.

Edit:  
I corrected my f-function, but unfortunately the solution is still a vector of NaNs

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [April 12, 2021, 4:40pm UTC](https://discourse.julialang.org/t/solving-system-of-nonlinear-equations-using-modelingtoolkit-jl-and-nonlinearsolve-jl/59096/5 "2021-04-12T16:40:02Z")

</div>

The issue is that with the choice of units, an order-1 guess for c[1] causes cosh and sinh to overflow. (At first I thought these must be CGS units, but in that case, `g = 9.81e02` cm/s^2, right?)

```julia
julia> m, g, L0 = 257*2.7, 9.81e03, 600e03;

julia> x1, x2, y1, y2 = 0, 500e03, 0, 0;

julia> c = [1;1;1];

julia> m*g/c[1]*(x1 + c[2])
6.807159000000001e6

julia> cosh(m*g/c[1]*(x1 + c[2]))
Inf

```

You can prevent the initial overflow by setting the guess for c1 to something large, e.g.

```julia
julia> c = [1e06;1;1];

julia> cosh(m*g/c[1]*(x1 + c[2]))
452.14957457816104

```

But a nonlinear solver is going to have a hard time dealing with a search space where there are large differences of scale on the components of c. Better to nondimensionalize the system or choose units where your numbers are order 1 or order 10. MKS should do.

---

<div class="post-metadata">

**Author:** ![wst](https://avatars.discourse-cdn.com/v4/letter/w/e79b87/32.png) [@wst](https://discourse.julialang.org/u/wst)\
**Post date:** [April 12, 2021, 6:10pm UTC](https://discourse.julialang.org/t/solving-system-of-nonlinear-equations-using-modelingtoolkit-jl-and-nonlinearsolve-jl/59096/6 "2021-04-12T18:10:10Z")

</div>

Thanks for pointing out the cosh-function shooting to infinty, guess I lack some intuition with hyperbolas.

The unit system I used was MMGS (realized I made an conversation error in my `m` term 😅).  
Switched to MKS for now, but as you said the nonlinear solver doesn’t like the conditioning (is this the right term here?) of the system and throws an `LinearAlgebra.SingularException(3)` error, even with `c_1`-guess of `10000`. I guess I’ll try to nondimensionalize the system next.

Edit:  
I played arround with my `c_1`-guess and at `c_1 = 10` I’m getting a somewhat realistic solution. Thanks again for your help @John_Gibson !

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [May 11, 2021, 11:53am UTC](https://discourse.julialang.org/t/solving-system-of-nonlinear-equations-using-modelingtoolkit-jl-and-nonlinearsolve-jl/59096/7 "2021-05-11T11:53:35Z")

</div>

If you rewrite the argument of `cosh` and `sinh` as follows

```julia
(m g / c1) (x1 + c2) = (L m g / c1) (x1/L + c2/L) 

                      = C1 (x1/L + C2)

```

with `C1 = Lmg/c1` and `C2 = c2/L`, then `C1` and `C2` should be order-1 constants. The nonlinear solve over `C1, C2` will be better-behaved and will avoid that overflow for order-1 initial guesses.

The right term here is preconditioning: transforming the problem so that it’s better suited to numerical solution methods.

---

<div class="post-metadata">

**Author:** ![wst](https://avatars.discourse-cdn.com/v4/letter/w/e79b87/32.png) [@wst](https://discourse.julialang.org/u/wst)\
**Post date:** [May 16, 2021, 10:14am UTC](https://discourse.julialang.org/t/solving-system-of-nonlinear-equations-using-modelingtoolkit-jl-and-nonlinearsolve-jl/59096/8 "2021-05-16T10:14:18Z")

</div>

Thanks for the clarification.  
I’ll try your rewritten arguments out as soon as I find some time for it.