python - numpy 数组更精确

Question

我有以下要评估的简单函数。

def f0(wt):
    term1 = (1 + np.cos(wt)**2) * (1 / 3 - 2 / (wt)**2)
    term2 = np.sin(wt)**2
    term3 = 4 / (wt)**3 * np.cos(wt) * np.sin(wt)
    return 0.5 * (term1 + term2 + term3)

对于 wt 的小值(1e-4 及以下的顺序)，我似乎在计算函数时遇到了数值问题。事实上，term1 和 term3 的值非常大且几乎相反，但是 term2 非常小。

我想我通过将 3 项的总和分成两部分来稍微改进一下，如下所示

def f1(wt):
    # Split the calculation to have more stability hopefully
    term1 = (1 + np.cos(wt)**2) * (1 / 3 - 2 / (wt)**2)
    term2 = np.sin(wt)**2
    term3 = 4 / (wt)**3 * np.cos(wt) * np.sin(wt)
    partial = term1 + term3
    return 0.5 * (partial + term2)

但是，对于非常小但正的wt值，我认为仍然存在数值问题。我希望此函数对于 wt 的任何正值都是平滑的，但是，正如您从附图中看到的那样，在低于 1e-3 的值时，存在异常伪像.

我的问题是:如果我已经在使用数据类型float64，我该如何提高 Numpy 的数值精度？

注意:我使用的是 64 位 Windows 10 计算机。我在其他 Stack Overflow 线程上读到类 np.float128 不可用。

完整代码片段

import numpy as np
import matplotlib.pyplot as plt

wt = np.logspace(-6, 1, 1000)

def f0(wt):
    term1 = (1 + np.cos(wt)**2) * (1 / 3 - 2 / (wt)**2)
    term2 = np.sin(wt)**2
    term3 = 4 / (wt)**3 * np.cos(wt) * np.sin(wt)
    return 0.5 * (term1 + term2 + term3)

def f1(wt):
     # Split the calculation to have more stability hopefully
     term1 = (1 + np.cos(wt)**2) * (1 / 3 - 2 / (wt)**2)
     term2 = np.sin(wt)**2
     term3 = 4 / (wt)**3 * np.cos(wt) * np.sin(wt)
     partial = term1 + term3
     return 0.5 * (partial + term2)

plt.figure()
plt.loglog(wt, f0(wt), label='f0')
plt.loglog(wt, f1(wt), label='f1')
plt.grid()
plt.legend()
plt.xlabel('wt')
plt.show()

Answer 1

用泰勒级数的前几项替换 sin 和 cosin 怎么样？然后 sympy 能够为您提供一个简单的结果，希望它更适合数值。

首先我稍微改变了你的函数，所以它给我一个 sympy 表达式。

from sympy import *
t = symbols('t')

def f0(wt):
    term1 = (1 + sympy.cos(wt)**2) * (sympy.Rational(1,3) - 2 / (wt)**2)
    term2 = sympy.sin(wt)**2
    term3 = 4 / (wt)**3 * sympy.cos(wt) * sympy.sin(wt)
    return sympy.Rational(1,2)*(term1 + term2 + term3)
expr = f0(t)
expr

现在我用泰勒多项式替换 sin 和 cos。

def taylor(f, n):
    return sum(t**i/factorial(i) * f(t).diff(t, i).subs(t,0) for i in range(n))

tsin = taylor(sin, 7)
tcos = taylor(cos, 7)

expr2 = simplify(expr.subs(sin(t),tsin).subs(cos(t),tcos))
f1 = lambdify(t, expr2, 'numpy')
expr2

最后我完全使用您的代码绘制它。请注意，我正在使用 sympys 选项制作一个 numpy ufunc。

wt = np.logspace(-6, 1, 1000)
plt.figure()
plt.loglog(wt, f0(wt), label='f0')
plt.loglog(wt, f1(wt), label='f1')
plt.grid()
plt.legend()
plt.xlabel('wt')
plt.show()

显然，此函数仅适用于零附近，对于 1 到 10 之间的值，您应该采用原始函数。但是，如果您需要说服力并且不关心替换泰勒多项式的函数看起来不错，您可以将度数提高到 25，使其在视觉上至少与您的函数一致，至少直到 10。

您可以组合这些函数，这样它就可以像这样用我的函数计算零附近的值，用你的函数计算另一个值。

def f2(wt):
    cond = np.abs(wt) > 1/10
    return np.piecewise(wt, [cond, ~cond], [f0,f1])