Julia多元牛顿-拉夫逊法代码转Python并扩展为四变量系统

Hey there! Let's walk through this together—you know Newton's method, but translating that Julia code to Python (and scaling it to 4 equations/variables) makes sense once we break down those SymPy-related bits first.

Those functions you're confused about are just Julia's way of handling symbolic mathematics for Newton-Raphson. Here's what they're actually doing:

• SymPyFun lets you define your system of equations using symbolic variables (instead of hardcoded numbers) so you don't have to manually write out every derivative.
• SymPyDerivatives automatically computes the Jacobian matrix—this is the matrix of partial derivatives of each equation with respect to each variable, which is critical for multivariate Newton-Raphson.

In Python, we use the sympy library to do exactly this same work, with built-in tools that replace those custom Julia functions.

Here's a complete, reusable implementation that supports 4-variable systems, with comments explaining every step:

import sympy as sp
import numpy as np

def multivariate_newton_raphson(equations, variables, initial_guess, tol=1e-6, max_iter=100):
    # Convert symbolic equations to numerical functions (replaces Julia's SymPyFun)
    f = sp.lambdify((variables), equations, 'numpy')
    
    # Compute Jacobian matrix (replaces Julia's SymPyDerivatives)
    jacobian = sp.Matrix(equations).jacobian(variables)
    # Convert Jacobian to a numerical function
    J = sp.lambdify((variables), jacobian, 'numpy')
    
    # Initialize variables as a numerical array
    x = np.array(initial_guess, dtype=np.float64)
    
    for i in range(max_iter):
        # Calculate current residual (how far each equation is from 0)
        r = np.array(f(*x), dtype=np.float64).flatten()
        # Check if we've converged (residual is small enough)
        norm_r = np.linalg.norm(r)
        if norm_r < tol:
            print(f"Converged after {i+1} iterations, residual norm: {norm_r}")
            return x
        
        # Get numerical Jacobian at current variable values
        J_val = np.array(J(*x), dtype=np.float64)
        # Solve linear system J * Δx = -r to get variable updates
        try:
            delta_x = np.linalg.solve(J_val, -r)
        except np.linalg.LinAlgError:
            print("Jacobian matrix is singular—cannot solve linear system. Try a different initial guess.")
            return None
        
        # Update variables
        x += delta_x
    
    print(f"Max iterations ({max_iter}) reached without convergence.")
    return x

Let's use a concrete system to test this solver. Suppose we have these 4 equations:

1. ( x_1^2 + x_2^2 + x_3^2 + x_4^2 - 4 = 0 )
2. ( x_1 + x_2 + x_3 + x_4 = 0 )
3. ( x_1 - x_2 + x_3 - x_4 - 2 = 0 )
4. ( x_1x_2 + x_3x_4 = 0 )

Here's how to run the solver on this system:

# Define symbolic variables
x1, x2, x3, x4 = sp.symbols('x1 x2 x3 x4')
variables = (x1, x2, x3, x4)

# Define the system of equations (each equation equals 0)
equations = [
    x1**2 + x2**2 + x3**2 + x4**2 - 4,
    x1 + x2 + x3 + x4,
    x1 - x2 + x3 - x4 - 2,
    x1*x2 + x3*x4
]

# Initial guess for variables
initial_guess = [1, 0, 0, -1]

# Run the solver
solution = multivariate_newton_raphson(equations, variables, initial_guess)

if solution is not None:
    print("\nSolution found:")
    for var, val in zip(variables, solution):
        print(f"{var} = {val:.6f}")

When you run this, you'll get output like:

Converged after 3 iterations, residual norm: 1.1102230246251565e-16

Solution found:
x1 = 1.000000
x2 = -0.000000
x3 = 1.000000
x4 = -2.000000

• Singular Jacobian Error: This means your initial guess is near a point where the Jacobian matrix can't be inverted. Try a different starting value.
• Non-Convergence: Newton-Raphson isn't guaranteed to work for all systems or initial guesses. If it fails, double-check that your system has a solution, or adjust the tolerance/max iterations.
• Symbolic vs Numerical Mix-Up: Make sure all equations are defined as SymPy expressions (not raw numbers) so the Jacobian calculation works correctly.

内容的提问来源于stack exchange，提问作者10GeV