如何用Python求解多解三角函数方程？寻求合适模块与算法

需要求解的三角函数方程为：
$$\tan(Nk) = \frac{d\sin k}{1 + d\cos k}$$
已知根据参数Δ的取值，方程在区间$(0, \pi)$内存在$M$个或$M-1$个解。此前尝试的采样+容差判断方法对容敏度极高，无法稳定工作，原代码如下：

import numpy as np

def f(k, N, d):
    return np.tan(k*(N)) - d*np.sin(k)/(1+d*np.cos(k))

k = np.linspace(0, np.pi, 10000+1, endpoint=False)[1:] 

def ksolutions(k,N,d):
    solutions=[]
    tol=4*1e-4
    for i in k:
        if abs(f(i,N,d)) < tol:
            solutions.append(i)
    
    print(solutions) 

可行的Python模块与算法 #

1. SciPy数值根查找模块 #

利用scipy.optimize.root_scalar针对每个连续区间单独查找根，是最可靠的方案：

• 原方程中$\tan(Nk)$在$k = \frac{(2m+1)\pi}{2N}$（$m=0,1,...,N-1$）处存在间断，因此可将$(0, \pi)$划分为$N$个连续子区间：$(0, \frac{\pi}{2N}), (\frac{\pi}{2N}, \frac{3\pi}{2N}), ..., (\frac{(2N-1)\pi}{2N}, \pi)$
• 在每个子区间内，函数连续且单调，适合用二分法（bisect）或牛顿法（newton）查找根

示例代码：

import numpy as np
from scipy.optimize import root_scalar

def f(k, N, d):
    return np.tan(N*k) - d*np.sin(k)/(1 + d*np.cos(k))

def find_solutions(N, d):
    solutions = []
    # 生成所有间断点，划分区间
    break_points = [(2*m + 1)*np.pi/(2*N) for m in range(N)]
    intervals = [(0, break_points[0])]
    for i in range(N-1):
        intervals.append((break_points[i], break_points[i+1]))
    intervals.append((break_points[-1], np.pi))
    
    for a, b in intervals:
        # 检查区间端点函数值是否接近零
        fa = f(a, N, d)
        fb = f(b, N, d)
        if abs(fa) < 1e-10:
            solutions.append(a)
            continue
        if abs(fb) < 1e-10:
            solutions.append(b)
            continue
        # 函数值异号则区间内存在根
        if np.sign(fa) != np.sign(fb):
            res = root_scalar(f, args=(N, d), method='bisect', bracket=[a, b], rtol=1e-10)
            if res.converged:
                solutions.append(res.root)
    return solutions

# 示例调用
N = 3
d = 0.5
solutions = find_solutions(N, d)
print("解：", solutions)

2. SymPy符号计算模块 #

如果需要解析解或高精度数值解，可使用SymPy：

• 用sympy.nsolve针对每个区间进行数值求解，或尝试sympy.solve获取解析解（但对于一般N和d可能无法得到闭合解）

示例代码：

import sympy as sp

k_sym = sp.symbols('k')
N = 3
d = 0.5

# 定义符号方程
eq = sp.tan(N*k_sym) - d*sp.sin(k_sym)/(1 + d*sp.cos(k_sym))

# 针对每个区间数值求解
solutions = []
break_points = [(2*m + 1)*sp.pi/(2*N) for m in range(N)]
intervals = [(0, break_points[0])]
for i in range(N-1):
    intervals.append((break_points[i], break_points[i+1]))
intervals.append((break_points[-1], sp.pi))

for a, b in intervals:
    try:
        sol = sp.nsolve(eq, k_sym, (a + b)/2)
        # 验证解是否有效
        if a < sol < b and abs(eq.subs(k_sym, sol)) < 1e-10:
            solutions.append(float(sol))
    except Exception:
        continue

print("解：", solutions)

原方法的问题 #

原代码通过均匀采样+容差判断的方式存在明显缺陷：

• 采样密度不足时会漏解，密度过高则效率低下
• 容差设置难以兼顾：容差太大容易误判非根点，容差太小则会漏掉接近根的点
• 无法处理函数在间断点附近的剧烈变化，容易出现错误判断

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