使用scipy.optimize.fsolve求解隐式方程时遇除零警告的解决方法

我在边坡稳定性分析中需要求解隐式方程的安全系数Fs，用自定义递归函数已得到结果1.8430，但调用scipy.optimize.fsolve时出现除零警告，推测是迭代过程中Fs取到0导致的。相关代码如下：

import numpy as np
from scipy.optimize import fsolve

# 定义方程函数
def fn4(Fs):
    return ( np.sum( (lamL * c + lamW * np.tan(np.radians(phi))) 
                 / (np.cos(lamA) * (1 + np.tan(lamA) * np.tan(np.radians(phi)) / Fs))
               ) / np.sum(lamW * np.sin(lamA))
           )

# 参数定义
c = 10
phi = 30
lamela1 = {'W':887.36, 'alpha':np.radians(46.71), 'L':19.325}
lamela2 = {'W':1624.8, 'alpha':np.radians(22.054), 'L':14.297}
lamela3 = {'W':737.43, 'alpha':np.radians(1.9096), 'L':13.258}
lamele = [lamela1, lamela2, lamela3]

lamW = np.array([[i['W'] for i in lamele]])
lamA = np.array([[i['alpha'] for i in lamele]])
lamL = np.array([[i['L'] for i in lamele]])

# 递归求解函数
def iterf(f, Fsi):
    if np.isclose(Fsi, fn4(Fsi)) == True:
        return Fsi
    else:
        return iterf(f, fn4(Fsi))

# 递归求解结果
print(iterf(fn4, 1))  # 返回1.8430

调用fsolve(fn4, 1)时出现警告：RuntimeWarning: divide by zero encountered in true_divide，原因是迭代过程中Fs可能被取到0或负数，触发分母为0的情况。以下是正确的解决方法：

1. 重构残差函数，避免Fs出现在分母 #

fsolve要求输入的是残差函数（即方程等于0的形式），而原函数fn4是迭代更新的表达式（Fs_new = fn4(Fs_old)），所以隐式方程应为Fs - fn4(Fs) = 0。更关键的是要把Fs从分母位置移出来，彻底避免除零：

def fn4_residual(Fs):
    # 整理方程，将Fs移到分子，消除分母中的Fs
    term1 = np.sum( (lamL * c + lamW * np.tan(np.radians(phi))) / np.cos(lamA) ) * Fs
    term2 = np.sum( lamW * np.sin(lamA) * (Fs + np.tan(lamA) * np.tan(np.radians(phi))) )
    return term1 - term2

# 调用fsolve，初始值设为1.5（更接近实际解，减少异常迭代）
result = fsolve(fn4_residual, 1.5)
print("fsolve求解结果:", result[0])  # 输出约1.8430

2. 给Fs添加最小值保护（保留原函数结构） #

如果不想大幅修改原函数，可以在计算时给Fs设置一个极小的下限，防止其变为0或负数：

def fn4_safe(Fs):
    # 限制Fs最小值为1e-6，避免除零
    Fs = max(Fs, 1e-6)
    return ( np.sum( (lamL * c + lamW * np.tan(np.radians(phi))) 
                 / (np.cos(lamA) * (1 + np.tan(lamA) * np.tan(np.radians(phi)) / Fs))
               ) / np.sum(lamW * np.sin(lamA))
           )

# 定义残差函数
def residual(Fs):
    return fn4_safe(Fs) - Fs

# 调用fsolve
result = fsolve(residual, 1)
print("fsolve求解结果:", result[0])

3. 选择更适合的求解器 #

可以使用scipy.optimize.root并指定method='lm'（勒文贝格-马夸尔特法），该方法对非线性方程的稳定性更好，尤其适合有物理约束的问题：

from scipy.optimize import root

sol = root(fn4_residual, 1.5, method='lm')
print("root求解结果:", sol.x[0])

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