Solve system of nonlinear equations w/ Newton-Raphson iteration

New to Julia
,
1

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)

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.
    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 w/ @StefanKarpinski @jw3126 @tkf.
    It seems to have fizzled out.
2 
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
3

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
4 
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
5

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
6

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

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

8

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
9

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

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.

11

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