Python数值方程根

工程分析中遇到的一个常见问题如下：给定函数$f(x)$，确定$f(x)=0$的$x$的值。解（$x$的值）被称为方程$f(x)=0$的根。

在继续进行之前，回顾一下函数的概念可能会有所帮助。 等式

$$y=f(x)$$

包含三个元素：输入值$x$，输出值$y$和计算$y$的规则$f$。 如果指定了规则$f$，则称指定该函数。 在数值计算中，规则始终是计算机算法。 它可能是一个函数语句，例如

$$f(x)=\cosh (x)\cos (x)-1$$

或包含数百或数千行代码的复杂过程。 只要算法为每个输入$x$生成输出$y$，它就可以作为函数

方程的根可以是实数，也可以是复数。 复数根很少计算，因为它们几乎没有物理意义。 多项式方程式是一个例外

$$a_0+a_{1}x+a_1{x}^2+\cdots +a_n{x}^n=0$$

复数根可能有意义的地方（例如，在阻尼振动分析中）。

通常，方程式可以具有任意数量的（实数）根或根本没有根。 例如，

$$\sin x-x=0$$

有一个根，即$x=0$，而

$$\tan x-x=0$$

有无限数量的根$(x=0,\pm 4.493,\pm 7.725,\dots )$。

所有找到根的方法都是需要起点（即，对根的估计）的迭代过程。 这个估计至关重要。 错误的起始值可能无法收敛，或者可能收敛到“错误”的根（与所求根不同的根）。 尚无用于估计根值的通用方法。 如果方程与物理问题相关联，则问题的上下文（物理洞察力）可能会建议根的大概位置。 否则，您可以对根进行系统的数值搜索。绘制函数是定位根的另一种方法，但这是无法编程的视觉过程。

增量搜索法

根的近似位置最好通过绘制函数来确定。 通常，基于一些点的非常粗糙的图就足以提供合理的起始值。 用于检测和括弧的另一个有用工具是增量搜索方法。 它也可以适用于计算根，但是这样做并不值得，因为本文中介绍的其他方法对于完成该任务更为有效。

增量搜索方法背后的基本思想很简单：如果$f\left(x_1\right)$和$f\left(x_2\right)$具有相反的符号，则间隔$\left(x_1,x_2\right)$中至少有一个根。 如果间隔足够小，则很可能包含单个根。 因此，可以通过在间隔$\Delta x$处评估函数并寻找符号变化来检测$f(x)$的零。

增量搜索方法存在一些潜在的问题：

• 如果搜索增量$\Delta x$大于两个根的间距，则有可能会丢失两个紧密间隔的根。
• 不会检测到双根（两个根重合）。
• $f(x)$的某些奇点（极点）可能被误认为是根。 例如，$f(x)=\tan x$会在$x=\pm \frac{1}{2}n\pi ,n=1,3,5,\dots$处改变符号，如图1所示。 但是，这些位置不是真正的零，因为该函数不与$x$轴交叉。

图1

算法代码

rootsearch函数在间隔$(a,b)$中以$dx$的增量搜索用户提供的函数$f(x)$的零。 如果搜索成功，则返回根的边界$(\mathrm{x}1,\mathrm{x}2)$；否则，返回根的边界。$\mathrm{x}1=\mathrm{x}2=\text{ None }$表示未检测到根。 在检测到第一个根（最接近$a$的根）之后，可以再次调用rootsearch并用$x2$代替，以查找下一个根。 只要rootsearch检测到根，就可以重复执行此操作。

## module rootsearch
’’’ x1,x2 = rootsearch(f,a,b,dx).
		Searches the interval (a,b) in increments dx for
		the bounds (x1,x2) of the smallest root of f(x).
		Returns x1 = x2 = None if no roots were detected.
’’’
from numpy import sign
def rootsearch(f,a,b,dx):
		x1 = a; f1 = f(a)
		x2 = a + dx; f2 = f(x2)
		while sign(f1) == sign(f2):
				if x1 >= b: return None,None
				x1 = x2; f1 = f2
				x2 = x1 + dx; f2 = f(x2)
		else:
				return x1,x2

示例1

$x^3-10x^2+5=0$的根位于间隔$(0,1)$中。 使用rootsearch函数可以四位数精度计算此根。

为了获得四位数的精度，我们需要不大于$\Delta x=0.0001$的搜索增量。 因此，以增量$x$搜索间隔（0，1）将需要10,000个功能评估。 下面的程序通过在四个阶段中接近根，将函数评估的数量减少到40个，每个阶段涉及10个搜索间隔（因此需要10个函数评估）。

二分法

在间隔$\left(x_1,x_2\right)$中将$f(x)=0$的根括弧后，可以使用几种方法对其进行接近。 对分方法通过将间隔连续减半直到间隔变得足够小来实现此目的。 此技术也称为间隔减半方法。 二等分不是用于计算根的最快方法，但它是最可靠的方法。 一旦将一个方括号括起来，二等分将始终在其上闭合。

二等分的方法使用与增量搜索相同的原理：如果区间$\left(x_1,x_2\right)$中有根，则$f\left(x_1\right)$和$f\left(x_2\right)$具有相反的符号。 为了将间隔减半，我们计算$f\left(x_3\right)$，其中$x_3=\frac{1}{2}\left(x_1+x_2\right)$是间隔的中点。 如果$f\left(x_2\right)$和$f\left(x_3\right)$具有相反的符号，则根必须在$\left(x_2,x_3\right)$中，并且我们通过用$x_3$替换原始绑定$x_1$来记录它。 否则，根位于$\left(x_1,x_3\right)$中，在这种情况下，$x_2$将替换为$x_3$。 无论哪种情况，新间隔$\left(x_1,x_2\right)$都是原始间隔的一半。 重复二等分，直到间隔减小到较小的值$\epsilon$，于是

$$\left|x_2-x_1\right|\le \epsilon$$

计算达到规定的$\epsilon$所需的平分数量很容易。 一个二等分后，原始间隔$\Delta x$减少为$\Delta x/2$，两个二等分后为$\Delta x/2^2$，$n$个二等分后为$\Delta x/2^n$。设定$\Delta x/2^n=\epsilon$，对n求解，我们得到

$$n=\frac{\ln (\Delta x/\epsilon )}{\ln 2},(1)$$

由于$n$必须为整数，因此使用$n$的上限（$n$的上限是大于$n$的最小整数）。

算法代码

bisection函数使用二分法计算已知位于区间$\left(x_1,x_2\right)$中的$f(x)=0$的根。 根据公式（1）可计算出将间隔减小到tol所需的二分数量n。通过将$\text{ switch }=1$，我们强制例程检查$f(x)$的幅度是否随着每个间隔减半而减小。 如果不是，则可能是错误的（可能“根”根本不是根，而是一个极点），并且返回root = None。 由于并非总是需要此功能，因此默认值为$\text{ switch }=0$。韩数error.err用来终止程序。

示例2

使用二分法找到$x^3-10x^2+5=0$的根，该根位于$(0,1)$到四位精度之间（此问题已在示例1中通过rootsearch解决了）。 该过程涉及多少功能评估？

代码如下：

请注意，我们将$\epsilon =0.0001$（二分中$\mathrm{t}\mathrm{o}\mathrm{l}=1.0\mathrm{e}-4$）设置为将精度限制为四个有效数字。 结果是

依据等式（1）

$$n=\frac{\ln (|\Delta x|/\epsilon )}{\ln 2}=\frac{\ln (1.0/0.0001}{\ln 2}=13.29$$

因此，二分法的for循环中的函数求值数为$\lceil 13.29\rceil =14$。在子例程的开头有另外2个求值，总共有16个函数求值。

示例3

用二分法求出区间（0，20）中$f(x)=x-\tan x$的所有零。 使用函数rootsearch和bisection。

请注意，tan x是奇异的，并且在x =π/ 2，3π/ 2，…处更改正负号。为防止二等分将这些点误认为是根，我们将switch设置为1。根与奇异点的接近度是另一个潜在的问题，可以通过在根搜索中使用小$\Delta x$来缓解。 选择$\Delta x=0.01$，我们得出以下结果代码：

详情参阅 - 亚图跨际