Python fsolve/root求解结果与MATLAB不符的优化方法咨询

Let's break down why your Python solution isn't matching MATLAB's results, fix the critical code mismatches, and then cover optimization tips to get consistent, accurate outputs.

First: Fix the Direct Code Errors #

Side-by-side comparison of your Python and MATLAB code reveals several key mismatches that are causing the result differences:

1. Depth Unit Conversion Mismatch #

In MATLAB, you convert depth to imperial units with depth = depth/0.3048 before passing it to the solver. In Python, you created dx but never used it—you’re still passing the original un-converted depth (2000) to fun(), which drastically changes the exponential terms in your equations.

2. Typo in the Second Equation #

Your Python code incorrectly uses x[1] twice in the second equation of fun(), while MATLAB correctly references x(1) and x(2). This typo completely breaks the second equation’s logic.

3. Missing dh Multiplier in First Equation #

MATLAB’s first equation includes g_e*cosd(incl)*dh, but Python’s F[0] only has g_e*math.cos(...)—you forgot to multiply by deltah (dh).

4. Loop Range Inconsistency #

MATLAB loops from n:-1:1 (covers all n elements), but Python’s range(n-1,0,-1) skips the first element (i=0), so you’re only solving 9 cases instead of 10.

5. Extra Cosine Term in Exponential #

MATLAB uses theta_2*depth in the exponential term, but Python incorrectly adds math.cos(math.radians(incl)) to that term. This was an unnecessary, incorrect addition that skews results.

Corrected Python Code #

Here’s the fixed version that aligns perfectly with your MATLAB logic:

import numpy as np
import math
from scipy import optimize

def fun(x, A, theta_1, theta_2, depth, dh, g_e, incl):
    F = np.zeros(2)
    # Fix: Add dh multiplier, remove extra cos(incl) in theta_2's exponential
    F[0] = x[0] * np.exp(theta_1 * depth) + x[1] * np.exp(theta_2 * depth) + g_e * math.cos(math.radians(incl)) * dh
    # Fix: Correct x[0] instead of x[1] for the first term in F[1]
    F[1] = ((1 + theta_1/A) * x[0] * np.exp(theta_1 * depth)) + ((1 + theta_2/A) * x[1] * np.exp(theta_2 * depth)) + (g_e * dh * math.cos(math.radians(incl)))
    return F

n = 10
original_depth = 2000
# Match MATLAB's unit conversion for depth
depth_converted = original_depth / 0.3048
t_1 = 0.001063321803317305
t_2 = -0.000485917956755497

# Recreate inclination array to match MATLAB's 1x100 shape
incl_1 = np.zeros(30)
incl_2 = np.linspace(0, 30, 30)
incl_3 = np.linspace(30, 80, 40)
incl_main = np.concatenate([incl_1, incl_2, incl_3])

g = 0.025
A = 0.0008948453688318264
# Match MATLAB's deltah calculation (uses original depth, not converted)
deltah = original_depth / n

mn = np.zeros((n, 2))
x0 = np.array([0, 0])
# Fix loop range to cover all n elements (0-indexed)
for i in range(n-1, -1, -1):
    mn[i] = optimize.fsolve(fun, x0, args=(A, t_1, t_2, depth_converted, deltah, g, incl_main[i]))

print(mn)

Optimization Tips for Better Solver Performance #

Once the code aligns with MATLAB’s logic, if you need to refine convergence or accuracy, try these steps:

1. Match MATLAB’s Solver Method #

MATLAB’s fsolve uses the Levenberg-Marquardt method by default for square systems. In Python, explicitly set this method in optimize.root for consistency:

from scipy.optimize import root
options = {'maxiter': 1000, 'xtol': 1e-12}
result = root(fun, x0, args=(A, t_1, t_2, depth_converted, deltah, g, incl_main[i]), method='lm', options=options)
mn[i] = result.x

2. Adjust Initial Guess #

Instead of using [0,0] for every iteration, use the previous iteration’s result as the initial guess—this helps the solver converge faster, especially since you’re looping backwards:

x0 = np.array([0, 0])
for i in range(n-1, -1, -1):
    mn[i] = optimize.fsolve(fun, x0, args=(A, t_1, t_2, depth_converted, deltah, g, incl_main[i]))
    x0 = mn[i]  # Update guess for next iteration

3. Provide the Jacobian #

Supplying the function’s Jacobian (gradient) helps the solver converge more accurately and quickly. Here’s how to define and use it:

def jac(x, A, theta_1, theta_2, depth, dh, g_e, incl):
    J = np.zeros((2,2))
    exp1 = np.exp(theta_1 * depth)
    exp2 = np.exp(theta_2 * depth)
    J[0,0] = exp1
    J[0,1] = exp2
    J[1,0] = (1 + theta_1/A) * exp1
    J[1,1] = (1 + theta_2/A) * exp2
    return J

# Use with fsolve
mn[i] = optimize.fsolve(fun, x0, args=(...), fprime=jac)

# Or with root
result = root(fun, x0, args=(...), jac=jac, method='lm')

4. Tune Solver Tolerances #

Adjust precision parameters to match MATLAB’s default or stricter settings:

# For fsolve
mn[i] = optimize.fsolve(fun, x0, args=(...), xtol=1e-10, maxfev=1000)

# For root
options = {'xtol': 1e-12, 'maxiter': 2000}

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