计算科学 - 寻找超越方程的前 N 个根 - 吾爱随笔录

寻找超越方程的前 N 个根

计算科学优化 Python 特殊功能寻根

2021-11-27 12:46:29

我需要找到第一个 $n$ 超越方程的根

F (k) = J_{m}^{'} (k r) Y_{m}^{'} (k) - J_{m}^{'} (k) Y_{m}^{'} (k r)

$\begin{equation} F(k) = J_m'(kr)Y_m'(k)-J'_m(k)Y'_m(kr) \end{equation}$

对于整数值 $m$ 和任何 $r \in [0,1)$ 在哪里 $J'$ 和 $Y'$ 是第一类和第二类贝塞尔函数的导数。这是环形圆柱体的特征函数展开中遇到的（标准？）问题。使用贝塞尔函数的递归性质，我们可以在没有导数的情况下重写这个方程，如下面的代码所示。

为了找到根，我fsolve在 python 中使用标准的根查找算法，这需要初始猜测。我已经收集到第一个根的一个不错的猜测将是 $m$ ，我们要为其求根的贝塞尔函数的阶数。从图中我还了解到，根的间距有些均匀，因此在找到第二个根后，我可以找到其他猜测。

根的间距是非线性的 $r$ . 我确信有一种强大的方法可以帮助进行根猜测，或者也许有更好的方法来做这件事。

from scipy.special import jn, yn 
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import argrelextrema


# Recurrence relations to relate derivative of bessels to normal bessel
Jp = lambda m,x : 0.5*(jn(m-1,x)-jn(m+1,x))
Yp = lambda m,x : 0.5*(yn(m-1,x)-yn(m+1,x))


# r is the ratio of inner/ outer radii and is [0,1)
F = lambda k,m,r : Jp(m,k*r)*Yp(m,k)-Jp(m,k)*Yp(m,k*r)

plt.figure()
k_array = np.linspace(0.1,80,5000)
m = 5

r = 0.9

F_array = F(k_array,m,r)



guesses = np.asarray([m,m/(1-r)])

roots = fsolve(F,guesses,args=(m,r,))
#iroots = argrelextrema(F_array**2,np.less)[0]
#print iroots




plt.plot(k_array,F_array)
#plt.plot(k_array[iroots],np.zeros(len(iroots)),"ro")
plt.plot(roots,np.zeros(len(roots)),"ro")

3个回答

这是一个单一维度上的解析函数的求根问题。一种标准技术是近似 $F(k)$ 作为多项式 $F(k) \approx c_0 + c_1 k + c_2 k^2 + \cdots$ 使用切比雪夫近似，并以半解析的方式计算多项式的根，例如通过建立伴生矩阵并计算其特征值。（鉴于 $F(k)$ 是解析的，它由多项式很好地逼近，并且 Chebyshev 多项式逼近的收敛性是指数的。）使用此过程找到的根可以使用牛顿法选择性地细化。

上述算法在 MATLAB 包 Chebfun 中实现。该算法的参考来自 John P. Boyd。

Boyd, John P. “通过切比雪夫展开和多项式求根计算实数区间上的零点。” SIAM 数值分析杂志40.5 (2002): 1666-1682。

我将用他建议的方法的 python 实现来补充@Richard Zhang 的答案 (+1)。MATLAB 包Chebfun已部分移植到python. 实际上有两个版本可用：chebpy和pychebfun.

这是根查找过程的实现pychebfun（方法与类似chebpy）

from scipy.special import jn, yn 
import pychebfun

# Recurrence relations to relate derivative of bessels to normal bessel
Jp = lambda m,x : 0.5*(jn(m-1,x)-jn(m+1,x))
Yp = lambda m,x : 0.5*(yn(m-1,x)-yn(m+1,x))


# r is the ratio of inner/ outer radii and is [0,1)
F = lambda k,m,r : Jp(m,k*r)*Yp(m,k)-Jp(m,k)*Yp(m,k*r)


m, r = 5, 0.9

# obtain the chebyshev approximation of the function for a specific domain 
f_cheb = pychebfun.Chebfun.from_function(lambda x: F(x, m, r), domain = (0.1,80))

# obtain the roots of the function in this domain
roots = f_cheb.roots()


#plot function and roots
k_array = np.linspace(0.1,80,5000)
plt.plot(k_array, F(k_array, m, r), label = '$F$')
plt.plot(f_cheb.roots(), F(f_cheb.roots(), m, r), 'o', label = 'roots')
plt.ylim(ymin=-.1, ymax = 0.1)
plt.legend()
plt.grid()
plt.xlabel('$k$');

作为答案的补充，我想说直接离散微分方程可能非常有效。我在上一个问题中使用了这种方法。我使用有限差分并使用特征求解器来找到根。除了按照牛顿法进行的扰动分析之外。

那里的建议也是使用切比雪夫近似。

其它你可能感兴趣的问题

上一篇密集系统的 Krylov 子空间方法下一篇计算科学竞赛。为什么没有？

寻找超越方程的前 N ​​个根

寻找超越方程的前 N 个根