Solve system of nonlinear equations w/ Newton-Raphson iteration

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:

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)


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.
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?
It seems to have fizzled out.
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.

1 Like

This almost works:

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

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

1 Like

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

1 Like

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

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.

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

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

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.

2 Likes

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

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.

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

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

1 Like