Conquer Nonlinear Systems with Python and SciPy#

What Is a Linear System?#

A linear system of equations has the form: [ \begin{cases} a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n = b_1 \ a_{21} x_1 + a_{22} x_2 + \cdots + a_{2n} x_n = b_2 \ \vdots \ a_{m1} x_1 + a_{m2} x_2 + \cdots + a_{mn} x_n = b_m \end{cases} ] Such algebraic systems can be efficiently solved using matrix operations, with solutions derived from linear algebra—for instance, by computing (\mathbf{x} = A^{-1}\mathbf{b}).

What Is a Nonlinear System?#

A general nonlinear system of (m) equations in (n) unknowns is: [ \begin{cases} f_1(x_1, \ldots, x_n) = 0 \ f_2(x_1, \ldots, x_n) = 0 \ \vdots \ f_m(x_1, \ldots, x_n) = 0 \end{cases} ] Here, each (f_i) can involve polynomial terms, exponentials, logarithms, trigonometric functions, or any complicated combination thereof. Closed-form solutions are rare; hence numerical approaches become essential. Techniques to solve such systems often rely heavily on iterative methods, starting with an initial guess and improving that guess until the functions are “close enough�?to zero.

Key Differences#

1. Uniqueness of solution: Linear systems typically have a unique solution if the coefficient matrix is non-singular. Nonlinear systems may have multiple solutions, no solutions, or an infinite set of solutions, depending on the behavior of the functions.
2. Sensitivity to initial guesses: Nonlinear solvers usually need a reasonable initial guess to converge; the path taken can lead to different roots for the same system if multiple solutions exist.
3. Numerical stability: Methods for nonlinear systems can sometimes diverge or get trapped in local minima, requiring robust strategies to detect and mitigate potential pitfalls.

Fixed-Point Iteration#

One straightforward (though not always robust) approach is fixed-point iteration. Suppose you want to solve (f(x) = 0). You rearrange it to (x = g(x)). Then you iterate: [ x^{(k+1)} = g\bigl(x^{(k)}\bigr). ] Convergence is not guaranteed in all cases, but the method can work well for simple scenarios or when a good transformation (g) is chosen.

Newton’s Method#

One of the most famous methods for solving (f(x)=0) for a single variable is Newton’s Method (also extended to multiple variables): [ x^{(k+1)} = x^{(k)} - \frac{f\bigl(x^{(k)}\bigr)}{f’\bigl(x^{(k)}\bigr)}. ] For multidimensional problems, you use the gradient (Jacobian) matrix to guide the iterative updates: [ \mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} - J^{-1}\bigl(\mathbf{x}^{(k)}\bigr) , \mathbf{f}\bigl(\mathbf{x}^{(k)}\bigr). ] Newton’s Method converges very quickly if you start near a solution and if the Jacobian is well-conditioned.

Other Methods: Secant and Quasi-Newton#

• Secant Method: Approximates the derivative numerically, avoiding symbolic differentiation.
• Quasi-Newton Methods (Broyden, DF-SANE, etc.): Build up an approximation to the Jacobian (or Hessian) progressively and can be highly effective for higher-dimensional problems.

Explanation#

1. We first define the single-variable function func(x).
2. We choose an initial guess x0 = 1.0.
3. We call fsolve(func, x0) and get an array with the solution.
4. We print and check the result; ideally, (f(\text{root})) is very close to zero.

Using bisect#

If we happen to know an interval in which the function changes sign, we can use the bounding approach bisect. For a single variable, this is often ideal because it guarantees convergence if f(a)*f(b) < 0.

1
from scipy.optimize import bisect
2

3
# Must find an interval [a, b] such that f(a)*f(b) < 0
4
a, b = 0, 2
5
root_bisect = bisect(func, a, b)
6
print("Bisection method root =", root_bisect)

What Did We Do?#

1. We grouped the two unknowns (x, y) using the list [x, y].
2. Our system function returns a list [eq1, eq2].
3. We fed an initial guess [1,1] into fsolve.
4. The solver returned the approximate solution, which should be (4,1) or (1,4). Actually, with the initial guess [1,1], fsolve should find (4,1).

Observing Multiple Solutions#

By nature, the system (xy=4) and (x+y=5) has two solutions: (4,1) and (1,4). If you try changing the initial guess, you might find the other solution. This demonstrates a key feature of nonlinear systems: the existence of multiple valid solutions.

Choosing an Initial Guess#

The success of methods like fsolve can hinge on a good initial guess. For physically oriented problems, domain knowledge and bounding conditions can guide your guess. Alternatively, random or systematic scanning can be used when no strong knowledge is available.

Handling Divergence and Local Minima#

• Divergence: If the solution is not converging, try adjusting the initial guess or switch to a different solver.
• Local Minima: In some contexts (e.g., optimization or solving residuals), your solver may get trapped in a local minimum. Techniques like scanning multiple initial guesses or using global optimization methods can help.

Scaling of Variables#

If your variables span vastly different magnitudes (e.g., some in the range of (10^{-9}) and others in the range of (10^{5})), it can impede numerical stability. Good practice is to scale variables to an approximately similar range, or let specialized solvers handle the scaling.

Convergence Tolerances#

You can adjust solver parameters such as xtol (tolerance in the solution) or ftol (tolerance in the function). Setting them too tight might cause slow convergence or non-convergence; setting them too loose might result in an imprecise solution.

SciPy’s root Function#

While fsolve is often enough, root offers a more generalized interface. You can specify method options such as "hybr", "lm", or "broyden1".

1
from scipy.optimize import root
2

3
def system(vars):
4
    x, y = vars
5
    return [x*y - 4, x + y - 5]
6

7
initial_guess = [1, 1]
8
res = root(system, initial_guess, method='lm')
9
print("Solution with 'lm':", res.x)

The 'lm' method uses the Levenberg-Marquardt approach, often quite robust for smaller systems. You could also try 'broyden1' or 'krylov' for bigger, sparser systems.

Large or Sparse Systems#

For extremely large systems (e.g., partial differential equations discretized into thousands of unknowns), you may need specialized solvers. SciPy has some limited coverage, but for advanced PDE problems, libraries like petsc4py, PyTrilinos, or specialized PDE solvers could be more efficient.

Jacobian and Derivatives#

Providing an analytical Jacobian can drastically improve performance, especially for Newton-based methods. SciPy allows you to pass a jac parameter to certain functions. If an analytical form is tough to derive, you might use sympy for symbolic differentiation or rely on finite differences.

Using least_squares#

SciPy’s least_squares is perfect for this:

1
import numpy as np
2
from scipy.optimize import least_squares
3

4
# Artificial data
5
x_data = np.array([0, 1, 2, 3, 4], dtype=float)
6
y_data = np.array([1.0, 2.7, 7.3, 19.4, 53.0], dtype=float)  # Suppose this is alpha=1, beta=1.0, plus noise
7

8
def residuals(params, x, y):
9
    alpha, beta = params
10
    # Return vector of residuals
11
    return y - alpha * np.exp(beta * x)
12

13
initial_guess = [1.0, 1.0]
14
result = least_squares(residuals, initial_guess, args=(x_data, y_data))
15
alpha_est, beta_est = result.x
16
print("Estimated alpha =", alpha_est)
17
print("Estimated beta =", beta_est)
18
print("Residuals norm =", result.cost)

In this snippet:

1. We define an array of (x)-values and corresponding measured (y)-values.
2. We define a residuals function that returns the difference between model and data.
3. By minimizing the sum of squared residuals, we effectively solve a system of partial derivatives set to zero.
4. The output result.x contains the best-fit (\alpha) and (\beta).

Understanding the Underlying Nonlinear System#

While least_squares is specialized for least-squares optimization, under the hood it’s akin to solving a system of equations: the gradient of the cost function with respect to the parameters is set to zero. This approach generalizes to more complicated models where you might have dozens or even hundreds of parameters.