A basic Newton method prototype.

newton_method.py

# -*- coding: utf-8 -*-

"""Newton method to solve three equations.
x^2 + y^2 = 800
y^2 + z^2 = 1300
x^2 + z^2 = 1300
[x = 20, y = 20, z = 30]
"""

from typing import Sequence, Callable
from numpy import array, ndarray, concatenate
from numpy.linalg import inv

no_inv = True
atom = 1e-10
inf = float('inf')
ConsFunc = Callable[[ndarray], float]


def partial(f: ConsFunc, x: ndarray, ind: int) -> float:
    x_new = x.copy()
    x_new[ind] += atom
    return (f(x_new) - f(x)) / atom


def gaussian_elimination(m: ndarray) -> ndarray:
    n = len(m)
    for i in range(0, n):
        # Search for maximum in this column
        max_el = abs(m[i, 0])
        max_row = i
        for k in range(i + 1, n):
            if abs(m[k, i]) > max_el:
                max_el = abs(m[k, i])
                max_row = k
        # Swap maximum row with current row (column by column)
        for k in range(i, n + 1):
            m[max_row, k], m[i, k] = m[i, k], m[max_row, k]
        # Make all rows below this one 0 in current column
        for k in range(i + 1, n):
            c = -m[k, i] / m[i, i]
            for j in range(i, n + 1):
                if i == j:
                    m[k, j] = 0
                else:
                    m[k, j] += c * m[i, j]
    # Backward substitution
    x = array([0 for _ in range(n)], dtype=float)
    for i in range(n - 1, -1, -1):
        x[i] = m[i, n] / m[i, i]
        for k in range(i - 1, -1, -1):
            m[k, n] -= m[k, i] * x[i]
    return x


def newton_method(x: ndarray, cons: Sequence[ConsFunc]) -> None:
    diff = array([inf, inf, inf], dtype=float)
    while abs(diff.sum()) > atom:
        jacobian = array([
            [partial(f, x, i) for i in range(len(cons))] for f in cons
        ], dtype=float)
        f_big = -array([f(x) for f in cons], dtype=float)
        if no_inv:
            # [J][x_new - x] = -[F]
            # Un-squeeze, 1D to 2D array: A -> A[None, :]
            jacobian = concatenate((jacobian, f_big[None, :].T), axis=1)
            diff = gaussian_elimination(jacobian)
        else:
            # [x_new] = [x] - [J]^-1[F]
            diff = inv(jacobian) @ f_big
        x += diff


def main():
    x_ = array([5, 5, 5], dtype=float)
    newton_method(x_, [
        lambda x: x[0] * x[0] + x[1] * x[1] - 800,
        lambda x: x[1] * x[1] + x[2] * x[2] - 1300,
        lambda x: x[0] * x[0] + x[2] * x[2] - 1300,
    ])
    print(x_)


if __name__ == '__main__':
    main()