# Solve system of nonlinear equations w/ Newton-Raphson iteration

**URL:** <https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383>\
**Category:** New to Julia\
**Tags:** nlsolve, newton-raphson\
**Created:** [August 6, 2020, 1:37am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383 "2020-08-06T01:37:06Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [August 6, 2020, 1:37am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/1 "2020-08-06T01:37:06Z")

</div>

Consider a boring system of 4 nonlinear equations in 4 variables:  
a=bcd  
a+b=cd  
a+b+c=d  
a+b+c+d=1  
Solution: (a,b,c,d)=(\frac{1}{42}, \frac{1}{7}, \frac{1}{3}, \frac{1}{2})

My code:

```julia
using Zygote;
NR(x0) = x0 - (fp(x0)^-1) * f(x0)
Nest(f,x0,n) = (n == 0) ? x0 : (f∘Nest)(f,x0,n-1)
# warmup 1
f(x)=x^2-2; fp(x)=f'(x); x0=1;
Nest.(NR,x0,0:10), √2
# warmup 2
f(x)=sin(x)-1; x0=1; fp(x)=f'(x);
Nest.(NR,x0,0:10), π/2
# 4 eqn system
f1(x)=x[1] - x[2]*x[3]*x[4] ;
f2(x)=x[1]+x[2] - x[3]*x[4] ;
f3(x)=x[1]+x[2]+x[3] - x[4] ;
f4(x)=x[1]+x[2]+x[3]+x[4] - 1 ;
f(x)=[f1(x);f2(x);f3(x);f4(x)]
fp(x)=Zygote.forward_jacobian(f,x)[2] |> transpose;
x0=[1;1;1;1];
Nest(NR,x0,100)

```

The answer is correct **BUT** :

1. I hate using the clunky-ugly syntax: `x[1]+x[2]+x[3]`  
How can I make it look like the math: `a+b+c`  
Is this a problem w/ the derivative package?
2. How can I get Julia to give me the answer in rational form?  
Is that even possible w/o symbolic math? I hope so.
3. I had to manually define the function `Nest(f,x0,n)` which is [routine in math software](https://reference.wolfram.com/language/ref/Nest.html).  
Function iteration (ie Dynamical Systems) show up in: Newton-Raphson, Mandelbrot set, value function iteration etc (all things used by Julians)  
Would it make sense for base Julia to include this?  
We had a [discussion about this on Github](https://github.com/JuliaLang/julia/pull/34251#issuecomment-635664141) w/ @StefanKarpinski @jw3126 @tkf.  
It seems to have fizzled out.

---

<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:** [August 6, 2020, 2:46am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/2 "2020-08-06T02:46:23Z")

</div>

```julia
julia> f1( (a, b, c, d) ) = a - b * c * d
f1 (generic function with 1 method)

julia> f1( [1, 2, 3, 4] )
-23

```

(Note the double parentheses in the definition of `f1`.)

An iterative method like Newton–Raphson converges (maybe) to the correct result, but never actually gets there, so I think it’s hopeless to expect the true rational result without symbolic methods.

The convention in Julia is that functions begin with lower-case letters, so `nest`.

`nest.(..., 0:10)` is inefficient since you’re recalculating everything each time.

There’s some way to write this with higher-order functions but I don’t remember how.

---

<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:** [August 6, 2020, 2:51am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/3 "2020-08-06T02:51:02Z")

</div>

This almost works:

```julia
julia> result = Nest(NR,x0,100)
4-element Array{Float64,1}:
 0.023809523809523815
 0.14285714285714288
 0.33333333333333337
 0.5

julia> rationalize.(result)
4-element Array{Rational{Int64},1}:
 62387527305566//2620276146833771
              1//7
              1//3
              1//2

```

---

<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:** [August 6, 2020, 2:55am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/4 "2020-08-06T02:55:01Z")

</div>

```julia
julia> result = Nest(NR, big.(x0), 100)
4-element Array{BigFloat,1}:
 0.02380952380952380952380952380952380952380952380952380952380952380952380952380974
 0.1428571428571428571428571428571428571428571428571428571428571428571428571428568
 0.3333333333333333333333333333333333333333333333333333333333333333333333333333348
 0.50

julia> rationalize.(ans)
4-element Array{Rational{Int64},1}:
 1//42
 1//7
 1//3
 1//2

```

Of course this does not really tell you that this rational solution is correct, but at least it gives a suggestion that you can then check to see if it really is a solution.

---

<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:** [August 6, 2020, 2:56am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/5 "2020-08-06T02:56:12Z")

</div>

But normally a system of nonlinear equations will not have a rational solution at all! This method will always report some rational “solution”.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [August 6, 2020, 3:07am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/6 "2020-08-06T03:07:46Z")

</div>

Even w/o `rationalize()` julia seems to know `d=.5` and not `.499999`

```julia
julia> Nest(NR,x0,100)
4-element Array{Float64,1}:
 0.023809523809523815
 0.14285714285714288
 0.33333333333333337
 0.5

```

You’re absolutely correct, the sol to most nonlinear systems will not be rational. In fact the sol to most diffeqs & integrals have no known closed forms. This is why I had to learn numerical methods in grad school.

My point is:

1. Numbers in Julia have a heierarchy of types: Real, Integer, Rational…
2. Suppose I want the “highest” number type possible in a solution, it would be nice if Julia can give me a solution (1.23333456, 1//2, 4) (Real, Rational, Integer) if possible.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [August 6, 2020, 3:16am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/7 "2020-08-06T03:16:45Z")

</div>

> [@dpsanders](#):
>
> `f1( (a, b, c, d) ) = a - b * c * d`

This would be awesome, but I’m not having luck differentiating this function of a tuple…

---

<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:** [August 6, 2020, 3:21am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/8 "2020-08-06T03:21:35Z")

</div>

Don’t know about Zygote, but ForwardDiff works fine:

```julia
julia> f1( (a, b, c, d) ) = a - b * c * d
f1 (generic function with 1 method)

julia> using ForwardDiff
ForwardDiff.
julia> ForwardDiff.gradient(f1, [1, 2, 3, 4])
4-element Array{Int64,1}:
   1
 -12
  -8
  -6

```

Julia can’t give you 1//2 and 4 automatically, since it has no way of knowing if the floating-point result of Newton-Raphson is “really” equal to those non-floating-point numbers or not.

In principle you could do Newton-Raphson with `Rational{BigInt}` all the way through (probably for fewer iterations). This will give you the exact result after N iterations, which will be a horrible enormous fraction which never becomes equal to 1//2, which it will converge to only in the limit of an infinite number of iterations.

In a more complicated problem Julia will return (for example) 0.4999999999999 instead of 0.5, and neither Julia nor you will have any idea whether it should “really” be 1//2 or not.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [August 6, 2020, 3:34am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/9 "2020-08-06T03:34:45Z")

</div>

@dpsanders thanks for the help.

1. I will try to rewrite w/ ForwardDiff
2. I’m not gonna worry about returning a desired type (int, rational) for now…
3. I’m still keen on making it easier to iterate functions in Julia.  
I [suggested](https://github.com/JuliaLang/julia/pull/34251#issuecomment-635664141) the Dynamical Systems notation f^n(x\_0) := nest(f,x\_0,n), bc it’s closest to the math & used A LOT in numerical math.  
At the very least Julia could use: `nest(f,x0,n), nestlist(f,x0,n)` and other routine [function iteration](https://reference.wolfram.com/language/guide/FunctionalIteration.html) methods.

---

<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:** [August 6, 2020, 3:43am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/10 "2020-08-06T03:43:32Z")

</div>

> [@Albert\_Zevelev](#):
>
> function iteration

I don’t recall hearing this called “nest” before; maybe that’s unique to Mathematica? It doesn’t bring to my mind iterated functions. I would just call it e.g. `iterated` or `orbit` or `trajectory`. Surely DynamicalSystems.jl has something like this implemented already.

---

<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:** [August 6, 2020, 5:31am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations-w-newton-raphson-iteration/44383/11 "2020-08-06T05:31:00Z")

</div>

> [@Albert\_Zevelev](#):
>
> I hate using the clunky-ugly syntax: `x[1]+x[2]+x[3]`  
> How can I make it look like the math: `a+b+c`  
> Is this a problem w/ the derivative package?

Use ModelingToolkit’s NonlinearSystem? [https://mtk.sciml.ai/dev/tutorials/nonlinear/](https://mtk.sciml.ai/dev/tutorials/nonlinear/)

> [@Albert\_Zevelev](#):
>
> I had to manually define the function `Nest(f,x0,n)` which is [routine in math software](https://reference.wolfram.com/language/ref/Nest.html).  
> Function iteration (ie Dynamical Systems) show up in: Newton-Raphson, Mandelbrot set, value function iteration etc (all things used by Julians)  
> Would it make sense for base Julia to include this?  
> We had a [discussion about this on Github](https://github.com/JuliaLang/julia/pull/34251#issuecomment-635664141) w/ @StefanKarpinski @jw3126 @tkf.  
> It seems to have fizzled out.

No, that’s generally not a great way to do it anyways.