使用Python在有限定义区间数值求解含参非线性方程组

Got it, let's tackle this problem step by step. You're looking to find a point (t₀, z₀) where two functions (one parameterized by t) match both in value and their z-derivative—with the key catch that at least one function is only defined over a specific interval. Since you already have a reasonable initial guess (tₛ, zₛ), here's how I'd approach this in Python:

First, frame your problem as a 2-dimensional nonlinear system of equations. We need two conditions to hold simultaneously:

1. ( F_1(t, z) = f(t, z, *args) - g(z, *args) = 0 )
2. ( F_2(t, z) = \frac{\partial f}{\partial z}(t, z, *args) - \frac{\partial g}{\partial z}(z, *args) = 0 )

Our goal is to find (t, z) that makes both ( F_1 ) and ( F_2 ) equal to 0, while respecting the domain constraints of your functions.

First, write a function that takes a combined array [t, z] and returns the two residuals ( F_1 ) and ( F_2 ). If you don't have analytical expressions for the derivatives, use numerical differentiation (Scipy has built-in tools for this).

Here's a concrete example with dummy functions (replace these with your actual functions):

import numpy as np
from scipy.optimize import least_squares
from scipy.misc import derivative

# Replace these with your actual f and g functions
def f(t, z, *args):
    return t * np.exp(-z) + args[0]  # Example: parameterized exponential

def g(z, *args):
    return np.sqrt(z) + args[1]  # Example: square root (only defined for z ≥ 0)

def objective(t_z, *args):
    t, z = t_z
    
    # Calculate first residual: f(t,z) - g(z)
    f_val = f(t, z, *args)
    g_val = g(z, *args)
    residual_1 = f_val - g_val
    
    # Calculate second residual: df/dz - dg/dz (numerical derivative)
    df_dz = derivative(lambda z_val: f(t, z_val, *args), z, dx=1e-6)
    dg_dz = derivative(lambda z_val: g(z_val, *args), z, dx=1e-6)
    residual_2 = df_dz - dg_dz
    
    return [residual_1, residual_2]

Since at least one function has a restricted domain, we need to keep the optimizer from wandering into undefined territory. There are two easy ways to handle this:

Option 1: Hard Bounds (for interval constraints) #

If your variables have clear upper/lower limits (e.g., ( z ≥ 0 ), ( t ∈ [0, 20] )), use the bounds parameter in least_squares:

def objective_with_bounds(t_z, *args):
    return objective(t_z, *args)

# Define bounds: ([t_min, z_min], [t_max, z_max])
bounds = ([0, 0], [20, np.inf])

Option 2: Penalty for Out-of-Domain Values #

If your constraints are more complex (e.g., a non-interval condition), add a large penalty to the residuals when variables are outside the valid domain:

def objective_with_penalty(t_z, *args, t_bounds=(0, 20), z_bounds=(0, np.inf)):
    t, z = t_z
    
    # Check if variables are in valid domain
    if not (t_bounds[0] <= t <= t_bounds[1]) or not (z_bounds[0] <= z <= z_bounds[1]):
        return [1e10, 1e10]  # Large penalty to push optimizer away
    
    return objective(t_z, *args)

Use your initial guess (tₛ, zₛ) to solve the system. least_squares is ideal here because it handles constraints well and is robust to noisy residuals:

# Example inputs
args = (2.0, 1.0)  # Extra arguments for f and g
initial_guess = [5.0, 2.0]  # Your (tₛ, zₛ) initial point

# Run the solver
result = least_squares(
    objective_with_penalty,
    initial_guess,
    args=args,
    bounds=bounds  # Use this if you went with hard bounds
)

# Check and print results
if result.success:
    t0, z0 = result.x
    print(f"Found solution: t₀ = {t0:.4f}, z₀ = {z0:.4f}")
    
    # Verify the solution
    f_val = f(t0, z0, *args)
    g_val = g(z0, *args)
    df_dz = derivative(lambda z: f(t0, z, *args), z0, dx=1e-6)
    dg_dz = derivative(lambda z: g(z, *args), z0, dx=1e-6)
    
    print(f"Validation: f - g = {f_val - g_val:.6f}, df/dz - dg/dz = {df_dz - dg_dz:.6f}")
else:
    print(f"Solver failed: {result.message}")

If you can write explicit expressions for ( \frac{\partial f}{\partial z} ) and ( \frac{\partial g}{\partial z} ), replace the numerical derivatives with analytical ones. This will boost accuracy and speed:

# Replace these with your actual analytical derivatives
def df_dz_analytic(t, z, *args):
    return -t * np.exp(-z)  # Derivative of our example f(t,z)

def dg_dz_analytic(z, *args):
    return 0.5 / np.sqrt(z)  # Derivative of our example g(z)

def objective_analytic(t_z, *args):
    t, z = t_z
    residual_1 = f(t, z, *args) - g(z, *args)
    residual_2 = df_dz_analytic(t, z, *args) - dg_dz_analytic(z, *args)
    return [residual_1, residual_2]

• Adjust the numerical derivative step size (dx) if you see unstable results—try values between ( 10^{-6} ) and ( 10^{-8} ).
• If the solver gets stuck, try a different initial guess or switch to a different method in scipy.optimize.root (like trust-exact or lm).
• Always validate the solution by plugging it back into your original functions and derivatives!

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