如何利用SciPy高效求解隐式方程的固定点迭代法？
- 内容介绍
- 相关推荐
本文共计743个文字，预计阅读时间需要3分钟。
使用`scipy.optimize.root_scalar`求解隐式方程如`(f_s=f(f_s))`的根时，可以采用以下步骤：
在边坡稳定性分析等工程计算中，常需求解形如 ( F_s = f(F_s) ) 的隐式方程——这类问题本质上是寻找函数的固定点（fixed point），而非零点（root）。初学者易误将固定点函数 fn4(Fs) 直接传入 scipy.optimize.fsolve，例如 fsolve(fn4, 1)，但这实际是在求解 fn4(Fs) = 0，不仅语义错误，更因 fsolve 内部可能尝试非物理域参数（如 Fs ≈ 0）而触发分母为零警告（RuntimeWarning: divide by zero），导致数值失败。
正确做法是：将原问题转化为标准的标量根查找问题，即定义残差函数 residual(Fs) = fn4(Fs) - Fs，再对其求零点。推荐使用 scipy.optimize.root_scalar ——它是专为单变量标量方程设计的现代接口，比已标记为“legacy”的 fsolve 更稳健、API 更清晰，并支持多种算法（如 'brentq'、'bisect'、'newton'）及显式搜索区间约束。
以下为完整可运行示例（已适配原始参数）：
import numpy as np from scipy.optimize import root_scalar # 定义几何与材料参数 c = 10.0 phi = 30.0 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]]) # shape: (1, 3) lamA = np.array([[i['alpha'] for i in lamele]]) # shape: (1, 3) lamL = np.array([[i['L'] for i in lamele]]) # shape: (1, 3) # 原始隐式函数 fn4(Fs) —— 返回待匹配的 Fs 值 def fn4(Fs): # 注意：Fs 是标量，参与广播运算；确保分母不为零（Fs > 0 物理合理） denominator = np.cos(lamA) * (1 + np.tan(lamA) * np.tan(np.radians(phi)) / Fs) numerator = lamL * c + lamW * np.tan(np.radians(phi)) return np.sum(numerator / denominator) / np.sum(lamW * np.sin(lamA)) # ✅ 正确构造残差函数：求解 fn4(Fs) - Fs == 0 residual = lambda Fs: fn4(Fs) - Fs # 方案1：使用初始猜测 x0（推荐 'secant' 或 'newton' 法） result = root_scalar(residual, x0=1.0, method='secant', xtol=1e-8) # 方案2：指定有界区间 bracket（更鲁棒，强制使用二分或 Brent 法） # result = root_scalar(residual, bracket=(0.5, 3.0), method='brentq', rtol=1e-10) print(f"Factor of Safety (Fs) = {result.root:.6f}") print(f"Converged? {result.converged}") print(f"Residual at solution: {residual(result.root):.2e}")
输出示例：
Factor of Safety (Fs) = 1.843012 Converged? True Residual at solution: 1.2e-14
综上，求解此类工程隐式方程的核心在于问题建模的准确性：明确是固定点问题 → 转化为残差为零的根查找 → 选用 root_scalar 并合理设置搜索域。这不仅规避了数值陷阱，更提升了结果的物理可信度与代码可维护性。
本文共计743个文字，预计阅读时间需要3分钟。
使用`scipy.optimize.root_scalar`求解隐式方程如`(f_s=f(f_s))`的根时，可以采用以下步骤：
在边坡稳定性分析等工程计算中，常需求解形如 ( F_s = f(F_s) ) 的隐式方程——这类问题本质上是寻找函数的固定点（fixed point），而非零点（root）。初学者易误将固定点函数 fn4(Fs) 直接传入 scipy.optimize.fsolve，例如 fsolve(fn4, 1)，但这实际是在求解 fn4(Fs) = 0，不仅语义错误，更因 fsolve 内部可能尝试非物理域参数（如 Fs ≈ 0）而触发分母为零警告（RuntimeWarning: divide by zero），导致数值失败。
正确做法是：将原问题转化为标准的标量根查找问题，即定义残差函数 residual(Fs) = fn4(Fs) - Fs，再对其求零点。推荐使用 scipy.optimize.root_scalar ——它是专为单变量标量方程设计的现代接口，比已标记为“legacy”的 fsolve 更稳健、API 更清晰，并支持多种算法（如 'brentq'、'bisect'、'newton'）及显式搜索区间约束。
以下为完整可运行示例（已适配原始参数）：
import numpy as np from scipy.optimize import root_scalar # 定义几何与材料参数 c = 10.0 phi = 30.0 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]]) # shape: (1, 3) lamA = np.array([[i['alpha'] for i in lamele]]) # shape: (1, 3) lamL = np.array([[i['L'] for i in lamele]]) # shape: (1, 3) # 原始隐式函数 fn4(Fs) —— 返回待匹配的 Fs 值 def fn4(Fs): # 注意：Fs 是标量，参与广播运算；确保分母不为零（Fs > 0 物理合理） denominator = np.cos(lamA) * (1 + np.tan(lamA) * np.tan(np.radians(phi)) / Fs) numerator = lamL * c + lamW * np.tan(np.radians(phi)) return np.sum(numerator / denominator) / np.sum(lamW * np.sin(lamA)) # ✅ 正确构造残差函数：求解 fn4(Fs) - Fs == 0 residual = lambda Fs: fn4(Fs) - Fs # 方案1：使用初始猜测 x0（推荐 'secant' 或 'newton' 法） result = root_scalar(residual, x0=1.0, method='secant', xtol=1e-8) # 方案2：指定有界区间 bracket（更鲁棒，强制使用二分或 Brent 法） # result = root_scalar(residual, bracket=(0.5, 3.0), method='brentq', rtol=1e-10) print(f"Factor of Safety (Fs) = {result.root:.6f}") print(f"Converged? {result.converged}") print(f"Residual at solution: {residual(result.root):.2e}")
输出示例：
Factor of Safety (Fs) = 1.843012 Converged? True Residual at solution: 1.2e-14
综上，求解此类工程隐式方程的核心在于问题建模的准确性：明确是固定点问题 → 转化为残差为零的根查找 → 选用 root_scalar 并合理设置搜索域。这不仅规避了数值陷阱，更提升了结果的物理可信度与代码可维护性。