5
$\begingroup$

Just as an exercise, I am numerically solving the following system of equations:

$$ \begin{equation} \begin{cases} x^2 + y^2 = 32 \\ 3x + 7y = 15 \end{cases} \end{equation} $$

I use the following code:

from scipy.optimize import root

nfev = 0
def fun(x):
    global nfev
    nfev += 1
    return x[0]**2 + x[1]**2 - 32, 3 * x[0] + 7 * x[1] - 15

njev = 0
def jac(x):
    global njev
    njev += 1
    jacobian = [[2 * x[0], 2 * x[1]], [3, 7]]
    return jacobian

x0 = [0, 0]
result = root(fun, x0, jac=jac)

print(result)
print(nfev, njev)

I get the following output:

 message: The solution converged.
 success: True
  status: 1
     fun: [-1.269e-08 -7.468e-10]
       x: [ 1.500e+00 -5.454e+00]
    nfev: 13
    njev: 1
    fjac: [[-9.147e-02 -9.958e-01]
           [-9.958e-01  9.147e-02]]
       r: [-9.918e+00  9.915e-01  1.082e+01]
     qtf: [-2.960e-08  4.167e-07]
15 2

There are 15 calls to fun and 2 calls to jac which are different from 13 and 1 for respectively nfev and njev result fields.

Removing jac=jac, I get the following output:

 message: The solution converged.
 success: True
  status: 1
     fun: [-3.553e-15  0.000e+00]
       x: [ 1.500e+00  5.454e+00]
    nfev: 17
    fjac: [[-7.626e-02 -9.971e-01]
           [ 9.971e-01 -7.626e-02]]
       r: [-1.003e+01 -8.319e-01  1.088e+01]
     qtf: [ 6.485e-10 -8.479e-09]
19 0

There are 19 calls to fun which is different than nfev = 17 shown in the result.

What could explain the discrepancies? Would providing the Jacobian accelerate the calculations or improve the accuracy? If so in which circumstances?

Adding two print statements at the end of the jac function as so:

print(x)
print(jacobian)

I get the following extra output:

[0 0]
[[0, 0], [3, 7]]
[0. 0.]
[[0.0, 0.0], [3, 7]]

This shows that jac is called twice with the same x values. This does not make much sense. Any clue?

$\endgroup$
3
  • $\begingroup$ This may also be a question specific to the default root solver, which is a modified Powell's method and is based on the implementation here math.utah.edu/software/minpack/minpack/hybrd.html $\endgroup$
    – whpowell96
    Nov 28, 2023 at 17:22
  • $\begingroup$ I tried to fix the math (the coefficients of the second equation don't agree with the code, which is really confusing), but "Edits must be at least 6 characters". sigh Hopefully somebody with god-like powers can fix the math. $\endgroup$
    – jeguyer
    Nov 28, 2023 at 19:20
  • 1
    $\begingroup$ @jeguyer Yeah, I fixed it. Thanks for the review $\endgroup$
    – Tarik
    Nov 29, 2023 at 6:07

2 Answers 2

9
$\begingroup$

At the point where you print out the Jacobian, adding

traceback.print_stack()

reveals that the first evaluation comes from _check_func(). This function validates the shape of the results of the user-supplied functions before calling the minpack nonlinear system solver.

$\endgroup$
2
$\begingroup$

Studying your code and comparing it with the provided system of nonlinear equations, there is an error in the following code snippet.

def fun(x):
    global nfev
    nfev += 1
    return x[0]**2 + x[1]**2 - 32, 3 * x[0] + 7 * x[0] - 15

Correction

def fun(x):
    global nfev
    nfev += 1
                                                 #x[0] - 15
    return x[0]**2 + x[1]**2 - 32, 3 * x[0] + 7 * x[1] - 15

Nonlinear system solution

With the correction made, your code generates the correct approximate solutions for the given problem. Since we have a nonlinear system formed by a line and a circle centered at the origin with radius $r = \sqrt{32}$, it admits two real solutions that satisfy the problem. For this problem, I wrote a code in JavaScript applying the Newton-Raphson algorithm. An important observation regarding this numerical method lies in the fact that we calculate the inverse of the Jacobian matrix; consequently, it cannot be null.

"use strict";

// This function performs the subtraction of two vectors.
function subtractionVector(A,B) {
    if (A.length !== B.length) {
        console.log('Problem! The inserted matrices are not the same size.');
        process.exit(0);
    } else {
        let matrixOut = [];
        for (let i = 0; i < A.length; i++) {
            matrixOut[i] = A[i] - B[i];
        }
        return matrixOut;
    }
}


// Multiplies an nxm matrix by a mx1 vector
function multiplyMatrix(A,B) {
    if (A[0].length !== B.length) {
        console.log("Problem! The inserted matrices have different orders.");
        process.exit(0);
    } else {
        let matrixOut = [];
        for (let i = 0; i < A.length; i++) {
            matrixOut[i] = 0;
            for (let j = 0; j < A[i].length; j++) {
                matrixOut[i] += A[i][j] * B[j];
            }
        }
        return matrixOut;
    }
}

// Inverse matrix
function inverseMatrix(A) {
    const n = A.length;
    for (let k = 0; k < n; k++) {
        const con = A[k][k];
        A[k][k] = 1;
        for (let j = 0; j < n; j++) {
            A[k][j] = A[k][j] / con;
        }
        for (let i = 0; i < n; i++) {
            if (i !== k) {
                const con = A[i][k];
                A[i][k] = 0.0;
                for (let j = 0; j < n; j++) {
                    A[i][j] = A[i][j] - A[k][j] * con;
                }
            }
        }
    }
    return A;
}

// Norm the vector
function norm(vector) {
    let sumSquare = 0;
    for (let i = 0; i < vector.length; i++) {
        sumSquare += vector[i]**2;
    }
    return Math.sqrt(sumSquare);
}

// Norm
function normInfinity(vector) {
    function valorMaximo(A){
        let n = A.length;
        let max = A[0];
        for (let i = 1; i<=n-1; i++) {
            if(max<A[i]){
                max = A[i];
            }
        }
        return max;
    }

    let matrixA = [];
    for (let i = 0; i < vector.length; i++) {
        matrixA[i] = Math.abs(vector[i]);
    }
    return valorMaximo(matrixA);
}

// System equation f(r) = [f1(x,y,...),f2(x,y,...)]
function f(r) {
    let x = r[0];
    let y = r[1];

    // System of equations
    let f1 = x**2+y**2-32;
    let f2 = 3*x+7*y-15;

    return [f1,f2];
}

// Jacobian system
// Attention! The Jacobian cannot be null, because in the
// iterative Newton-Rapshon method we calculate the inverse
// of this, which is a square matrix.
function Jacobian(r) {
    let x = r[0];
    let y = r[1];

    // Partial derivative of f1(x,y) with respect to x and y
    let f1_x = 2*x; // Partial derivative of x
    let f1_y = 2*y; // Partial derivative of y

    // Partial derivative of f2(x,y) with respect to x and y
    let f2_x = 3; // Partial derivative of x
    let f2_y = 7; // Partial derivative of y

    return [[f1_x,f1_y],[f2_x,f2_y]];
}

function NewtonRapshon(r0,tol) {
    let counter = 0, x = r0, s = r0;
    do {
        s = x; // previous
        x = subtractionVector(x,multiplyMatrix(inverseMatrix(Jacobian(x)),f(x)));
        if(counter>1000){
            console.log('It didn\'t converge!');
            break;
        }
        //console.log(normInfinity(subtractionVector(x,s)))
        counter++;
    } while (norm(subtractionVector(x,s))>tol);
    for (let i = 0; i < x.length; i++) {
        x[i] = parseFloat(x[i].toFixed(3));
    }
    console.log(`\nNumber of iterations: ${counter}`);
    return x;
}

// Tolerance
let tol = 1.0e-15; // Tolerance

// First root
let r0_1 = [4,5]; // First kick
let sol_1 = NewtonRapshon(r0_1,tol);
console.log('Solution: ',sol_1);

// Second root
let r0_2 = [-5,-3]; // First kick
let sol_2 = NewtonRapshon(r0_2,tol);
console.log('Solution: ',sol_2);

Final result

enter image description here

$\endgroup$
3
  • $\begingroup$ Thanks. I applied your correction. BTW, where do you run your javascript? Using node.js? Are there numerical libraries in JavaScript. I mean mainstream libraries. $\endgroup$
    – Tarik
    Nov 29, 2023 at 8:16
  • $\begingroup$ Hello, you don't need to download and install anything on your computer! Just use your preferred browser. Copy the JavaScript code, right-click on the open web page (any page), click 'Inspect.' When you click 'Inspect,' it will open a page called DevTools; click on 'Console,' paste the JavaScript code, and press enter. $\endgroup$ Nov 29, 2023 at 8:24
  • $\begingroup$ In terms of numerical libraries for numerical-scientific computation with JavaScript: Yes, they exist! However, Python, MATLAB, Octave, Scilab, and Julia are more well-known and widely used. I am familiar with one called Math.js (mathjs.org), which is quite similar to SciPy, but SciPy is much larger and more widely employed (at least that's what it seems). I only implemented it in JavaScript to apply the logic of the Newton-Raphson numerical method. $\endgroup$ Nov 29, 2023 at 8:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.