如何用scipy.fsolve限制变量范围求解纳米光纤HE11模式方程

我来帮你解决这个问题——你需要限制bsk的范围在n2到n1之间，因为你的目标函数在这个区间外要么无定义，要么会产生复数，导致fsolve无法正常工作。咱们一步步来调整你的代码：

先修正原代码的几个关键问题 #

首先你的代码里有几个容易导致错误的细节：

• 波长单位转换错误：纳米转米应该是lamb_m = lamb_nm * 1e-9，而不是lamb_nm*1，否则波数k的计算完全不对
• 变量名冲突：func2里的k0=k和你导入的修正贝塞尔函数k0重名，会覆盖函数引用，导致后续计算出错
• 初始猜测值可能踩边界：比如循环里初始bskg=1.000刚好等于n2，容易触发平方根的边界问题

方法一：修改目标函数，强制限制bsk范围 #

我们可以在目标函数里加入判断：如果bsk超出[n2, n1]区间，就返回一个极大的数值，让fsolve自动往有效区间内搜索。同时确保初始猜测值落在有效区间内。

修改后的完整代码：

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from scipy.special import jv, kn, jvp, kvp, j0, j1, k0, k1
from scipy.optimize import fsolve
import numpy as np
from numpy.lib import scimath as SM
import matplotlib.pyplot as plt

knp = kvp
kv = kn

def getBeta(a=125, n1=1.4, n2=1, lamb_nm=500, bskg=1.2):
    # 修正波长单位转换：纳米转米
    lamb_m = lamb_nm * 1e-9
    k = 2 * np.pi / lamb_m

    def func2(bsk):
        nu = 1
        # 重命名变量，避免和导入的k0函数冲突
        k0_val = k
        # 先判断bsk是否在有效区间内
        if bsk <= n2 or bsk >= n1:
            # 超出范围返回极大值，引导求解器往有效区间走
            return 1e10
        
        X = a * np.sqrt(k0_val**2 * (n1**2 - bsk**2))
        Y = a * np.sqrt((bsk**2 - n2**2) * k0_val**2)
        V = np.sqrt(X**2 + Y**2)
        
        # 计算目标函数
        term1 = -(n1**2 + n2**2)/(2*n1**2) * kvp(1, Y)/(Y * k1(Y))
        term2 = 1 / X**2
        sqrt_term = np.sqrt(
            ((n1**2 - n2**2)/(2*n1**2) * kvp(1, Y)/(Y * k1(Y)))**2 
            + (bsk**2 / n1**2) * (1/X**2 + 1/Y**2) * 2
        )
        term3 = j0(X)/(X * j1(X))
        
        result = term1 + term2 - sqrt_term - term3
        return result

    # 确保初始猜测值在有效区间内
    if not (n2 < bskg < n1):
        bskg = (n1 + n2) / 2  # 如果初始值无效，用区间中点替代
    beta = fsolve(func2, bskg)
    return beta[0]  # fsolve返回数组，取第一个元素

plt.figure()
plt.title("β vs Fiber Radius (HE11 Mode)")
x = np.arange(100, 1000)  # 光纤半径，单位nm
bskg = 1.1  # 初始猜测值落在1到1.4之间
y = np.array([], dtype=float)

for i in x:
    beta = getBeta(lamb_nm=500, a=i, bskg=bskg)
    y = np.append(y, beta)
    bskg = beta  # 用前一次的结果作为下一次的初始猜测，加速收敛

plt.plot(x, y)
plt.xlabel("Fiber Radius (nm)")
plt.ylabel("β (Propagation Constant)")
plt.show()

方法二：使用带边界约束的优化方法 #

fsolve本身不支持边界约束，如果你想更严格地限制bsk的范围，可以用scipy.optimize.minimize，把问题转化为最小化目标函数的平方，同时设置边界条件：

修改getBeta函数中的求解部分：

from scipy.optimize import minimize

def getBeta(a=125, n1=1.4, n2=1, lamb_nm=500, bskg=1.2):
    # ... 前面的代码和方法一一致，省略 ...
    
    # 定义最小化的目标：函数值的平方
    def min_func(bsk):
        return func2(bsk)**2
    
    # 设置边界：bsk必须在(n2, n1)之间，加小偏移避免边界上的复数问题
    bounds = [(n2 + 1e-6, n1 - 1e-6)]
    # 使用L-BFGS-B方法，支持边界约束
    res = minimize(min_func, x0=bskg, bounds=bounds, method='L-BFGS-B')
    
    if res.success:
        return res.x[0]
    else:
        # 如果求解失败，返回区间中点或者抛出警告
        print(f"求解失败: {res.message}")
        return (n1 + n2)/2

这种方法更可靠，因为求解器会严格在你指定的边界内搜索，不会出现超出范围的情况。

关键说明 #

• 有效区间：因为n1=1.4，n2=1，所以bsk必须满足1 < bsk < 1.4，否则X或Y会出现复数，导致目标函数无意义
• 初始猜测值：尽量落在有效区间内，这样求解器收敛更快，也更容易找到正确的解
• 波长单位：一定要确保单位转换正确，否则波数k的计算会完全错误，导致整个求解结果无效

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