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

Hey there,
I’m fairly new to Julia and this forum (in fact my first post here).
I’m currently trying to solve a system of 3 nonlinear equations that emerged as the solution of the hanging chain problem (following this 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:

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

1 Like

This is a @YingboMa question.

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> 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
``````
1 Like

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

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> 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> 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.

3 Likes

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 !

1 Like

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

``````(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.

2 Likes

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