# Techniques to improve resolution of a system of nonlinear equations

**URL:** <https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712>\
**Category:** Numerics\
**Created:** [October 20, 2020, 10:09pm UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712 "2020-10-20T22:09:39Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [October 20, 2020, 10:09pm UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/1 "2020-10-20T22:09:39Z")

</div>

Hello,

I would like to solve a set of N = 300 independent 1D nonlinear equations of the type f\_i(x\_i) = a\_i where the a\_i are given and each function f\_i is a univariate strictly increasing bijection. There is a unique solution for each nonlinear equation. So far, I am using NLsolve.jl with the Newton method to solve simultaneously the equations. I am already using the fact that the Jacobian is diagonal (all the problems are independent).

However the f\_i have some bumps (see Figure) that are sometimes challenging for the Newton method if the initial condition is far from the root. How can I improve the convergence?  
Is there a way to exploit that the Jacobian has only strictly positive values?

![monotone](https://global.discourse-cdn.com/julialang/original/3X/b/0/b01d215aa792cc4419d50ff51d895dfe9b07002d.png)

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [October 20, 2020, 10:57pm UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/2 "2020-10-20T22:57:54Z")

</div>

If the problems are independent, why don’t you just solve each equation separately?

---

<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:** [October 20, 2020, 11:17pm UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/3 "2020-10-20T23:17:37Z")

</div>

> [@mleprovost](#):
>
> Is there a way to exploit that the Jacobian has only strictly positive values?

Yes, since it’s 1D and monotonically increasing you can use a bracketing method like bisection if you find one value to the left and one value to the right.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 21, 2020, 12:21am UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/4 "2020-10-21T00:21:10Z")

</div>

Doesn’t Newton’s method converge faster than bisection (at least once you’re close)?

---

<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:** [October 21, 2020, 12:47am UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/5 "2020-10-21T00:47:48Z")

</div>

Bisection is much more stable though. The best thing to do would be to use a Falsi method like in [https://github.com/JuliaComputing/NonlinearSolve.jl](https://github.com/JuliaComputing/NonlinearSolve.jl) which is a mixture of bracketing for stability but derivative-based for acceleration.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [October 21, 2020, 12:54am UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/7 "2020-10-21T00:54:37Z")

</div>

The main reason is that the Newton method is much faster to converge than the bisection method.

I have added a line search and that seems to solve the problem.  
Thank you all for your recommendations.

I will check NonlinearSolve.jl

---

<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:** [October 21, 2020, 12:56am UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/8 "2020-10-21T00:56:34Z")

</div>

for reference if you’re not familiar with the method: [Regula falsi - Wikipedia](https://en.wikipedia.org/wiki/Regula_falsi#The_regula_falsi_(false_position)_method)

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [October 21, 2020, 1:50pm UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/9 "2020-10-21T13:50:11Z")

</div>

Bisection is much more robust, and is derivative-free.

For problems with “bumps” like this, plain vanilla Newton’s method can easily diverge. There are safeguards against this, but why bother when you have bisection.

Newton’s method (in its original, simple form) is not a practical _general_ method for univariate or multivariate problems. It is useful when you can prove reasonable convergence analytically, which happens for globally “nice” problems, which are quite rare.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [October 21, 2020, 2:51pm UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/10 "2020-10-21T14:51:19Z")

</div>

I thought the standard (?) method for 1d functions was Brent’s method, which combines bisection and Newton to give a method that is guaranteed to converge (if I remember correctly). I believe this is implemented in the Roots.jl package.

If you need guarantees then there’s IntervalRootFinding.jl

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [October 21, 2020, 3:10pm UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/11 "2020-10-21T15:10:30Z")

</div>

> [@dpsanders](#):
>
> Brent’s method, which combines bisection and Newton

Not quite — it [uses the secant method](https://en.wikipedia.org/wiki/Brent%27s_method), so it does not need derivatives directly.

Brent’s method can be better than bisection _for some problems_, but I think that first we should establish whether the problem requires a multivariate solver. If the problem is separable, the major efficiency gain will be from that, and the choice of univariate method is of secondary importance compared to it.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [October 21, 2020, 3:25pm UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/12 "2020-10-21T15:25:59Z")

</div>

Yes, the problem is separable (all the equations are independent), but my code works best if I perform the evaluation of all the f\_i(x\_i) at once. Also, the computation of the Jacobian is less expensive than the evaluation of the function.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [October 21, 2020, 4:56pm UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/13 "2020-10-21T16:56:06Z")

</div>

Start with bisection and turn on Newton when things start looking good (ie all the |f|s are small). This is the algorithm HP used years ago for scalar equations in their programmable calculators.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [October 22, 2020, 8:52am UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/14 "2020-10-22T08:52:44Z")

</div>

There is no technical reason that should prevent you from implementing a coordinate-wise bisection or other method, it is just a matter of bookkeeping, the algorithm is the same.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [October 22, 2020, 10:28am UTC](https://discourse.julialang.org/t/techniques-to-improve-resolution-of-a-system-of-nonlinear-equations/48712/15 "2020-10-22T10:28:36Z")

</div>

This is wise advice. You can use a scalar bisection + Newton-Armijo solver (or roll you own) and loop over the functions. In that way you are likely to converge for most of them and will be able to identify any corner cases for special attention. The line search is really important for problems like this.

One problem with using Newton-Armijo for systems to solve this problem is that the hardest of the equations will govern the line search for all of them. Take this example, please

f1(x) = atan(x); f2(x) = atan(10_x); f3(x) = atan(100_x);

with x0=1 as the initial iterate for all three equations.

f3 = 0 is a much harder problem for Newton-Armijo that f1=0, which is easy. Solving as a system will make all three problems appear to be hard.