数值评估¶
基础¶
可以使用 .evalf() 方法或 N() 函数将精确的 SymPy 表达式转换为浮点近似值（十进制数）。N(expr, <args>) 等同于 sympify(expr).evalf(<args>)。
>>> from sympy import *
>>> N(sqrt(2)*pi)
4.44288293815837
>>> (sqrt(2)*pi).evalf()
4.44288293815837
默认情况下，数值计算的精度为15位小数。你可以选择性地将所需的精度（应为正整数）作为参数传递给 evalf 或 N：
>>> N(sqrt(2)*pi, 5)
4.4429
>>> N(sqrt(2)*pi, 50)
4.4428829381583662470158809900606936986146216893757
支持复数：
>>> N(1/(pi + I), 20)
0.28902548222223624241 - 0.091999668350375232456*I
如果表达式包含符号或由于其他原因无法进行数值计算，调用 .evalf() 或 N() 将返回原始表达式，或在某些情况下返回部分计算的表达式。例如，当表达式是展开形式的多项式时，系数会被计算：
>>> x = Symbol('x')
>>> (pi*x**2 + x/3).evalf()
3.14159265358979*x**2 + 0.333333333333333*x
你也可以使用标准的 Python 函数 float()、complex() 将 SymPy 表达式转换为常规的 Python 数字：
>>> float(pi)
3.141592653589793
>>> complex(pi+E*I)
(3.141592653589793+2.718281828459045j)
如果使用这些函数，未能将表达式计算为显式数字（例如，如果表达式包含符号）将引发异常。
本质上没有上精度限制。例如，以下命令计算π/e的前100,000位：
>>> N(pi/E, 100000)
...
这显示了π的第999,951到1,000,000位数字：
>>> str(N(pi, 10**6))[-50:]
'95678796130331164628399634646042209010610577945815'
高精度计算可能会很慢。建议（但完全可选）安装 gmpy (https://github.com/aleaxit/gmpy)，这将显著加快上述计算的速度。
浮点数¶
SymPy 中的浮点数是 Float 类的实例。Float 可以通过第二个参数自定义精度来创建：
>>> Float(0.1)
0.100000000000000
>>> Float(0.1, 10)
0.1000000000
>>> Float(0.125, 30)
0.125000000000000000000000000000
>>> Float(0.1, 30)
0.100000000000000005551115123126
如最后一个示例所示，某些 Python 浮点数作为输入时仅精确到大约 15 位数字，而其他浮点数（分母为 2 的幂的那些，如 0.125 = 1/8）则是精确的。要从高精度十进制数创建 Float，最好传递一个字符串、Rational 或对 Rational 进行 evalf：
>>> Float('0.1', 30)
0.100000000000000000000000000000
>>> Float(Rational(1, 10), 30)
0.100000000000000000000000000000
>>> Rational(1, 10).evalf(30)
0.100000000000000000000000000000
一个数字的精度决定了 1) 在进行算术运算时使用的精度，以及 2) 打印数字时显示的位数。当两个精度不同的数字一起用于算术运算时，结果使用较高的精度。0.1 +/- 0.001 和 3.1415 +/- 0.0001 的乘积的不确定性约为 0.003，但显示的精度为 5 位。
>>> Float(0.1, 3)*Float(3.1415, 5)
0.31417
因此，显示的精度不应被用作误差传播或有效数字运算的模型；相反，这种方案被用来确保数值算法的稳定性。
N 和 evalf 可以用来改变现有浮点数的精度：
>>> N(3.5)
3.50000000000000
>>> N(3.5, 5)
3.5000
>>> N(3.5, 30)
3.50000000000000000000000000000
准确性和错误处理¶
当输入到 N 或 evalf 的是一个复杂的表达式时，数值误差传播成为一个问题。例如，考虑第100个斐波那契数和优秀的（但不精确的）近似值 \(\varphi^{100} / \sqrt{5}\)，其中 \(\varphi\) 是黄金比例。使用普通的浮点数算术，从这些数中减去对方会导致错误的完全抵消：
>>> a, b = GoldenRatio**1000/sqrt(5), fibonacci(1000)
>>> float(a)
4.34665576869e+208
>>> float(b)
4.34665576869e+208
>>> float(a) - float(b)
0.0
N 和 evalf 会跟踪错误并自动增加内部使用的精度，以便获得正确的结果：
>>> N(fibonacci(100) - GoldenRatio**100/sqrt(5))
-5.64613129282185e-22
不幸的是，数值评估无法区分一个表达式是否恰好为零，还是仅仅非常小。因此，工作精度默认被限制在大约100位数字。如果我们尝试计算第1000个斐波那契数，会发生以下情况：
>>> N(fibonacci(1000) - (GoldenRatio)**1000/sqrt(5))
0.e+85
返回的数字中缺少位数表明 N 未能达到完全的精度。结果表明表达式的量级小于 10^84，但这并不是一个特别好的答案。要强制使用更高的工作精度，可以使用 maxn 关键字参数：
>>> N(fibonacci(1000) - (GoldenRatio)**1000/sqrt(5), maxn=500)
-4.60123853010113e-210
通常，maxn 可以设置得非常高（数千位数），但请注意，在极端情况下这可能会导致显著的减速。或者，可以将 strict=True 选项设置为强制抛出异常，而不是静默地返回一个精度低于请求的值：
>>> N(fibonacci(1000) - (GoldenRatio)**1000/sqrt(5), strict=True)
Traceback (most recent call last):
...
PrecisionExhausted: Failed to distinguish the expression:
-sqrt(5)*GoldenRatio**1000/5 + 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
from zero. Try simplifying the input, using chop=True, or providing a higher maxn for evalf
如果我们添加一个术语，使得斐波那契近似变得精确（Binet公式的完整形式），我们会得到一个精确为零的表达式，但 N 并不知道这一点：
>>> f = fibonacci(100) - (GoldenRatio**100 - (GoldenRatio-1)**100)/sqrt(5)
>>> N(f)
0.e-104
>>> N(f, maxn=1000)
0.e-1336
在已知会发生此类取消的情况下，chop 选项很有用。这基本上会将一个数字的实部或虚部中的非常小的数值替换为精确的零：
>>> N(f, chop=True)
0
>>> N(3 + I*f, chop=True)
3.00000000000000
在需要去除无意义数字的情况下，重新评估或使用 round 方法是有用的：
>>> Float('.1', '')*Float('.12345', '')
0.012297
>>> ans = _
>>> N(ans, 1)
0.01
>>> ans.round(2)
0.01
如果你处理的是一个不包含浮点数的数值表达式，它可以被计算到任意精度。为了将结果相对于给定的十进制数进行四舍五入，round 方法是有用的：
>>> v = 10*pi + cos(1)
>>> N(v)
31.9562288417661
>>> v.round(3)
31.956
求和与积分¶
和（特别是无穷级数）与积分可以像常规的闭式表达式一样使用，并且支持任意精度的计算：
>>> var('n x')
(n, x)
>>> Sum(1/n**n, (n, 1, oo)).evalf()
1.29128599706266
>>> Integral(x**(-x), (x, 0, 1)).evalf()
1.29128599706266
>>> Sum(1/n**n, (n, 1, oo)).evalf(50)
1.2912859970626635404072825905956005414986193682745
>>> Integral(x**(-x), (x, 0, 1)).evalf(50)
1.2912859970626635404072825905956005414986193682745
>>> (Integral(exp(-x**2), (x, -oo, oo)) ** 2).evalf(30)
3.14159265358979323846264338328
默认情况下，使用 tanh-sinh 求积算法来评估积分。该算法对于光滑的被积函数（甚至是对端点奇异的积分）非常高效且稳健，但对于高度振荡或有中点不连续性的积分可能表现不佳。在许多情况下，evalf/N 将正确估计误差。对于以下积分，结果是准确的，但仅精确到四位小数：
>>> f = abs(sin(x))
>>> Integral(abs(sin(x)), (x, 0, 4)).evalf()
2.346
将这个积分分成两部分会更好：
>>> (Integral(f, (x, 0, pi)) + Integral(f, (x, pi, 4))).evalf()
2.34635637913639
一个类似的例子是以下振荡积分：
>>> Integral(sin(x)/x**2, (x, 1, oo)).evalf(maxn=20)
0.5
通过告诉 evalf 或 N 使用振荡积分算法，可以更有效地处理它：
>>> Integral(sin(x)/x**2, (x, 1, oo)).evalf(quad='osc')
0.504067061906928
>>> Integral(sin(x)/x**2, (x, 1, oo)).evalf(20, quad='osc')
0.50406706190692837199
振荡积分需要一个包含因子 cos(ax+b) 或 sin(ax+b) 的被积函数。注意，许多其他振荡积分可以通过变量变换转换为这种形式：
>>> init_printing(use_unicode=False)
>>> intgrl = Integral(sin(1/x), (x, 0, 1)).transform(x, 1/x)
>>> intgrl
oo
/
|
| sin(x)
| ------ dx
| 2
| x
|
/
1
>>> N(intgrl, quad='osc')
0.504067061906928
如果级数收敛得足够快，无穷级数使用直接求和。否则，使用外推法（通常是欧拉-麦克劳林公式，但也包括理查森外推法）来加速收敛。这使得可以对缓慢收敛的级数进行高精度评估：
>>> var('k')
k
>>> Sum(1/k**2, (k, 1, oo)).evalf()
1.64493406684823
>>> zeta(2).evalf()
1.64493406684823
>>> Sum(1/k-log(1+1/k), (k, 1, oo)).evalf()
0.577215664901533
>>> Sum(1/k-log(1+1/k), (k, 1, oo)).evalf(50)
0.57721566490153286060651209008240243104215933593992
>>> EulerGamma.evalf(50)
0.57721566490153286060651209008240243104215933593992
欧拉-麦克劳林公式也用于有限级数，使得它们可以在不计算所有项的情况下快速近似：
>>> Sum(1/k, (k, 10000000, 20000000)).evalf()
0.693147255559946
请注意，evalf 做了一些假设，这些假设并不总是最优的。为了对数值求和进行精细控制，手动使用 Sum.euler_maclaurin 方法可能是值得的。
对于有理超几何级数（其中项是多项式、幂、阶乘、二项式系数等的乘积），使用了特殊的优化。N/evalf 可以非常快速地将此类型的级数求和到高精度。例如，这个拉马努金关于 pi 的公式可以在几分之一秒内通过一个简单的命令求和到 10,000 位数：
>>> f = factorial
>>> n = Symbol('n', integer=True)
>>> R = 9801/sqrt(8)/Sum(f(4*n)*(1103+26390*n)/f(n)**4/396**(4*n),
... (n, 0, oo))
>>> N(R, 10000)
3.141592653589793238462643383279502884197169399375105820974944592307816406286208
99862803482534211706798214808651328230664709384460955058223172535940812848111745
02841027019385211055596446229489549303819644288109756659334461284756482337867831
...
数值简化¶
函数 nsimplify 尝试找到一个在数值上等于给定输入的公式。此功能可用于猜测近似浮点输入的精确公式，或猜测复杂符号输入的更简单公式。nsimplify 使用的算法能够识别简单分数、简单代数表达式、给定常数的线性组合，以及任何前述内容的某些基本函数变换。
可选地，nsimplify 可以传递一个包含常量的列表（例如 pi）和一个最小数值容差。以下是一些基本示例：
>>> nsimplify(0.1)
1/10
>>> nsimplify(6.28, [pi], tolerance=0.01)
2*pi
>>> nsimplify(pi, tolerance=0.01)
22/7
>>> nsimplify(pi, tolerance=0.001)
355
---
113
>>> nsimplify(0.33333, tolerance=1e-4)
1/3
>>> nsimplify(2.0**(1/3.), tolerance=0.001)
635
---
504
>>> nsimplify(2.0**(1/3.), tolerance=0.001, full=True)
3 ___
\/ 2
以下是几个更高级的示例：
>>> nsimplify(Float('0.130198866629986772369127970337',30), [pi, E])
1
----------
5*pi
---- + 2*e
7
>>> nsimplify(cos(atan('1/3')))
____
3*\/ 10
--------
10
>>> nsimplify(4/(1+sqrt(5)), [GoldenRatio])
-2 + 2*GoldenRatio
>>> nsimplify(2 + exp(2*atan('1/4')*I))
49 8*I
-- + ---
17 17
>>> nsimplify((1/(exp(3*pi*I/5)+1)))
___________
/ ___
1 / \/ 5 1
- - I* / ----- + -
2 \/ 10 4
>>> nsimplify(I**I, [pi])
-pi
----
2
e
>>> n = Symbol('n')
>>> nsimplify(Sum(1/n**2, (n, 1, oo)), [pi])
2
pi
---
6
>>> nsimplify(gamma('1/4')*gamma('3/4'), [pi])
___
\/ 2 *pi