数值计算#
基础#
精确的SymPy表达式可以转换为浮点近似值（十进制数），使用 .evalf() 方法或 N() 功能。 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.1415926535...
>>> complex(pi+E*I)
(3.1415926535...+2.7182818284...j)
如果使用这些函数，未能将表达式求值为显式数字（例如，表达式包含符号）将引发异常。
基本上没有精度上限。例如，以下命令计算π/e的前100000位：
>>> N(pi/E, 100000)
...
显示的是pi的999951到1000000位：
>>> str(N(pi, 10**6))[-50:]
'95678796130331164628399634646042209010610577945815'
High-precision calculations can be slow. It is recommended (but entirely optional) to install gmpy (https://github.com/aleaxit/gmpy), which will significantly speed up computations such as the one above.
浮点数#
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 或 evalf 一 Rational ：
>>> 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个Fibonacci数和极好的（但不是精确的）近似值 \(\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
如果处理的是不包含浮点的数值表达式，则可以将其计算为任意精度。要将结果相对于给定的小数舍入，舍入方法很有用：
>>> 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, wrap_line=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
如果无穷级数收敛速度足够快，则使用直接求和法。另外，外推方法（通常是Euler-Maclaurin公式和Richardson外推法）可以加快收敛速度。这允许对缓慢收敛的序列进行高精度评估：
>>> 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
Euler-Maclaurin公式也适用于有限级数，允许快速近似，而无需计算所有项：
>>> Sum(1/k, (k, 10000000, 20000000)).evalf()
0.693147255559946
注意 evalf 做出一些并不总是最优的假设。对于对数值求和的微调控制，手动使用该方法可能是值得的 Sum.euler_maclaurin .
有理超几何级数（其中的项是多项式、幂、阶乘、二项式系数等的乘积）使用了特殊的优化。 N/evalf 这种类型的求和级数非常快，精度很高。例如，使用一个简单的命令，可以在不到一秒的时间内将pi的这个Ramanujan公式相加为10000位数字：
>>> 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