数值积分方法quad导数错误及多参数优化适用性咨询

Let's break down your problems and work through solutions step by step:

1. Fixing the ValueError in Numerical Integration #

The error you're seeing comes from using scipy.integrate.quadrature—this method relies on computing derivatives of your integrand to construct orthogonal polynomials for integration, which fails because your complex symbolic function doesn't play nicely with its derivative calculation.

Solution: Switch to scipy.integrate.quad #

quad is an adaptive quadrature method that doesn't require derivative calculations, making it far more robust for complex numerical functions like yours. Additionally, you'll want to convert your SymPy symbolic function to a fast numerical lambda function using lambdify to avoid slow symbol-to-numeric conversions during integration.

Here's the modified code:

First, update your function definitions to generate a numerical integrand:

import numpy as np
import math as m
import scipy.integrate as spi
import sympy as sp

# Define constants
gammaee = 5.55e-6
MJpsi = 3.096916
alphaem = 1/137
lambdasq = 0.09
Ca = 3
qOsq = 2

def qbarsq(qsq): return (qsq + MJpsi**2)/4
def xx(qbarsq_val, w): return 4*qbarsq_val/(4*qbarsq_val - MJpsi**2 + w**2)

# Define symbolic variables
x, NN_sym, a_sym, b_sym, ktsq_sym, qbarsq_sym, w_sym = sp.symbols('x NN a b ktsq qbarsq w')

def xg(a, b, NN, ktsq, x): 
    return NN*(x**(-a))*(ktsq**b)*sp.exp(sp.sqrt((16*Ca/9)*sp.log(1/x)*sp.log((sp.log(ktsq/lambdasq))/(sp.log(qOsq/lambdasq)))))

def func(NN, a, b, x, ktsq): 
    return (-x*sp.diff(sp.log(xg(a, b, NN, ktsq, x)), x))

def Rg(NN, a, b, ktsq, x): 
    return 2**(2*func(NN,a,b,x,ktsq)+3)/sp.sqrt(sp.pi)*sp.gamma(func(NN,a,b,x,ktsq)+5/2)/sp.gamma(func(NN,a,b,x,ktsq)+4)

def FktsqDeriv(NN, a, b, x, ktsq): 
    return sp.diff(Rg(NN,a,b,ktsq,x)*xg(a,b,NN,ktsq,x), ktsq)

def Fktsq1(qbarsq, ktsq, NN, a, b, w): 
    return FktsqDeriv(NN,a,b,x,ktsq).subs(x, 4*qbarsq/(4*qbarsq - MJpsi**2 + w**2))

def fA_sym(qbarsq, ktsq, NN, a, b, w): 
    return Fktsq1(qbarsq,ktsq,NN,a,b,w) * (1/qbarsq) * (1/(qbarsq + ktsq))

# Convert symbolic function to fast numerical lambda
fA_num = sp.lambdify(
    (ktsq_sym, NN_sym, a_sym, b_sym, w_sym, qbarsq_sym), 
    fA_sym(qbarsq_sym, ktsq_sym, NN_sym, a_sym, b_sym, w_sym),
    modules=['numpy', 'math']
)

Then rewrite your integration function using quad:

def integrated_f(NN, a, b, w, qbarsq):
    lower = 1
    upper = (w**2 - MJpsi**2)/4
    # Ensure valid integration interval (lower < upper)
    if upper < lower:
        lower, upper = upper, lower
    # Use quad instead of quadrature
    return spi.quad(fA_num, lower, upper, args=(NN, a, b, w, qbarsq))

# Test the integration
a=0.1
NN=0.5
b=-0.2
w=89
qbarsq=5
result = integrated_f(NN,a,b,w,qbarsq)
print(f"Integral result: {result[0]}, Estimated error: {result[1]}")

This should eliminate the derivative-related ValueError, as quad doesn't rely on derivative calculations.

2. Feasibility of Using scipy.optimize.minimize with Nelder-Mead #

Yes, using Nelder-Mead for your parameter fitting is absolutely feasible—and actually a great choice! Here's why, plus key considerations:

Why Nelder-Mead works here #

Nelder-Mead is a gradient-free optimization method, which is perfect for your case because:

• Your objective function (the integral result) doesn't have an easy-to-compute analytical gradient.
• It handles non-smooth or complex functions well, as long as the function is continuous over the parameter space.

Key steps to implement it #

You'll need to wrap your integrated_f into an objective function that accepts a single parameter array (for NN, a, b) and returns a scalar value (the integral result, or a loss function if you're fitting to data).

Example code:

import scipy.optimize as spo

def objective(params, w, qbarsq):
    NN, a, b = params
    integral_val, _ = integrated_f(NN, a, b, w, qbarsq)
    # Replace this with your actual loss function if fitting to data
    # e.g., return np.abs(integral_val - experimental_data)
    return integral_val

# Initial parameter guess
x0 = [0.5, 0.1, -0.2]
# Fixed parameters for integration
w_fixed = 89
qbarsq_fixed = 5

# Run optimization
opt_result = spo.minimize(
    objective, 
    x0, 
    args=(w_fixed, qbarsq_fixed), 
    method='Nelder-Mead', 
    options={'max_iter': 1000}
)

print("Optimization results:")
print(opt_result)

Important notes #

1. Objective function scalar output: minimize requires your objective function to return a single scalar. If you're fitting to experimental data, compute a loss (e.g., squared error between your integral and data) instead of returning the integral directly.
2. Integration stability: Ensure your integration interval is always valid (lower < upper) for all parameter values during optimization—invalid intervals can cause NaNs or errors that break the optimization.
3. Speed: Numerical integration is computationally expensive, and Nelder-Mead requires many function evaluations. For faster results, consider:
   • Using a more efficient integration method if possible.
   • Adding bounds to your parameters (if you know valid ranges) to reduce the search space.
   • Trying other gradient-free methods like Powell if Nelder-Mead is too slow.

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