含非解析积分的超越方程数值求解代码报错排查

需要为数组C = np.linspace(0, 1000, 100)中的每个值，求解积分方程的解y：

∫₀ʸ [ (n(i)/(n(0)cosA))² -1 ]^(-0.5) di = C，其中n(i)=k*i+n₀

编写的Python代码运行时出现RuntimeWarning和IntegrationWarning，无法正常得到结果。

原代码如下：

import numpy as np
import scipy.integrate as integrate
from scipy.optimize import root
import matplotlib.pyplot as plt

n_0 = 1
k = 5
A = np.deg2rad(30)

def n(y):
    return k*y + n_0

def f(y):
    f = lambda i: (((n(i)) / (n(0) * np.cos(A)))**2 - 1)**(-0.5)
    result, error = integrate.quad(f, 0, y)
    return result

C = np.linspace(0, 1000, 1000)
dy = np.array([])

for arg in C: 
    def integral(y):
        return f(y) - arg
    dy = np.append(dy, root(integral, 1).x)

plt.plot(C, dy)
plt.show()

报错信息：

d:\FMF\MAF2\simulacija.py:14: RuntimeWarning: invalid value encountered in double_scalars
  f = lambda i: (((n(i)) / (n(0) * np.cos(A)))**2 - 1)**(-0.5)
d:\FMF\MAF2\simulacija.py:15: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  result, error = integrate.quad(f, 0, y)
d:\FMF\MAF2\simulacija.py:15: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  result, error = integrate.quad(f, 0, y)

核心问题分析 #

1. 被积函数定义域问题：求解器root在迭代过程中可能尝试负的y值，导致积分区间包含i<0的情况，此时n(i)=k*i+n0可能小于n0*cosA，使得被积函数中根号内为负数，触发invalid value警告。
2. 数值积分效率与稳定性问题：原代码使用数值积分求解，当积分区间较大时，quad会遇到细分次数不足、舍入误差等问题，触发IntegrationWarning，且数值积分的计算效率较低。
3. 循环内定义函数的冗余：每次循环都定义新的integral函数，既影响效率也可能引发变量作用域问题。

解决方案 #

步骤1：推导积分的解析解 #

原积分可以通过变量替换得到解析表达式，完全替代数值积分，避免积分相关的警告和效率问题：
推导过程：
令 ( u = \frac{k i + n_0}{n_0 \cos A} )，则 ( di = \frac{n_0 \cos A}{k} du )，积分上下限从 ( u=\frac{1}{\cos A} ) 变为 ( u=\frac{k y + n_0}{n_0 \cos A} )。
利用积分公式 ( \int \frac{1}{\sqrt{u^2-1}} du = \ln(u+\sqrt{u^2-1}) + C )，最终得到积分结果：
[
f(y) = \frac{n_0 \cos A}{k} \left[ \ln\left( \frac{k y + n_0 + \sqrt{(k y + n_0)^2 - (n_0 \cos A)^2}}{n_0 (1+\sin A)} \right) \right]
]

步骤2：优化方程求解逻辑 #

• 使用单变量求解器root_scalar替代root，更适合一维方程求解，且能指定搜索区间限制y>0，避免出现无效的负数值。
• 提前定义目标函数，避免循环内重复定义函数。

import numpy as np
from scipy.optimize import root_scalar
import matplotlib.pyplot as plt

n_0 = 1
k = 5
A = np.deg2rad(30)
cosA = np.cos(A)
sinA = np.sin(A)

# 积分的解析表达式
def integral_value(y):
    numerator = k * y + n_0
    sqrt_term = np.sqrt(numerator**2 - (n_0 * cosA)**2)
    log_arg = (numerator + sqrt_term) / (n_0 * (1 + sinA))
    return (n_0 * cosA / k) * np.log(log_arg)

# 针对每个C值，定义目标方程：integral_value(y) - C = 0
def target_func(y, C_val):
    return integral_value(y) - C_val

# 生成C数组
C = np.linspace(0, 1000, 1000)
dy = np.zeros_like(C)

# 遍历求解每个C对应的y
for idx, c_val in enumerate(C):
    # 使用root_scalar，指定搜索区间(0, 1e6)，避免负数
    sol = root_scalar(target_func, args=(c_val,), bracket=[1e-6, 1e6], method='brentq')
    dy[idx] = sol.root

# 绘图
plt.plot(C, dy)
plt.xlabel('C')
plt.ylabel('y')
plt.title('Solution y vs C')
plt.show()

• 解析解替代数值积分：直接使用推导后的解析公式计算积分值，彻底解决数值积分的警告和效率问题。
• 限制搜索区间：root_scalar的bracket参数限定y在正数区间，避免求解器尝试无效的负数值，消除invalid value警告。
• 高效求解：单变量求解器root_scalar比root更适合一维问题，计算更快更稳定。

内容的提问来源于stack exchange，提问作者Edward Henry Brenner