SymPy 数值计算完全指南evalf、N 与任意精度求值详解【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympySymPy 作为一套纯 Python 实现的计算机代数系统其数值求值Numerical Evaluation能力是连接精确符号世界与浮点计算世界的桥梁。本文以官方文档 doc/src/modules/evalf.rst 为骨架结合 sympy/core/evalf.py 等源码实现系统讲解evalf()方法与N()函数的使用、Float精度模型、误差处理机制、级数与积分的高精度数值求值以及nsimplify数值反推公式等实战技能。读完本文你将能熟练地对任意 SymPy 表达式进行从 15 位到上万位的高精度数值计算并正确应对精度耗尽、振荡积分、慢收敛级数等棘手场景。一、求值基础evalf()与N()1.1 两种等价入口精确的 SymPy 表达式可以通过.evalf()方法或N()函数转换为浮点近似值。两者的关系在 sympy/core/evalf.py 的N函数定义中一目了然def N(x, n15, **options): return sympify(x, rationalTrue).evalf(n, **options)即N(expr, args)等价于sympify(expr).evalf(args)区别仅在于N会先用rationalTrue对输入做符号化处理这能保证输入字符串中的浮点数在求值时以精确有理数参与运算。同时.n()方法与.evalf()也是完全等价的别名见源码中n evalf一行。基本用法 from sympy import * N(sqrt(2)*pi) 4.44288293815837 (sqrt(2)*pi).evalf() 4.442882938158371.2 指定精度位数默认情况下数值求值精确到 15 位十进制数字。可以通过第二个参数传入期望的精度正整数 N(sqrt(2)*pi, 5) 4.4429 N(sqrt(2)*pi, 50) 4.4428829381583662470158809900606936986146216893757从源码 EvalfMixin.evalf 可以看到精度n会先经dps_to_prec(n)换算为二进制精度prec内部求值则在prec 4的保护位下进行最后再舍入到用户要求的位数。这也是N(expr, 50)尾数完全可信的原因。1.3 复数与符号表达式复数是天然支持的 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这一行为的实现细节在 evalf 内部函数evalf()首先在evalf_table中按表达式类型查找专用的求值器Add、Mul、Pow、exp、sin、Integral、Sum等均有注册见 _create_evalf_table查不到时回退到_eval_evalf当子结果仍不是数值时抛出的NotImplementedError会被上层捕获并返回原表达式见 evalf 方法 的异常回退逻辑。1.4 与 Python 原生数值互转也可以使用标准 Python 函数float()、complex()将 SymPy 表达式转换为普通 Python 数值 float(pi) # doctest: SKIP 3.141592653589793 complex(piE*I) # doctest: SKIP (3.1415926535897932.718281828459045j)注意使用这些函数时若表达式无法求值为显式数值例如含有符号会抛出异常——这与.evalf()的失败时返回原式策略截然不同。1.5 几乎没有上限的精度SymPy 的数值求值在原理上不设精度上限。下面的命令计算 π/e 的前 100,000 位数字 N(pi/E, 100000) #doctest: SKIP ...下面的命令展示 π 小数点后第 999,951 到 1,000,000 位 str(N(pi, 10**6))[-50:] #doctest: SKIP 95678796130331164628399634646042209010610577945815高精度计算可能较慢。官方文档建议可选安装 gmpy 以显著加速此类计算——SymPy 通过 sympy/external/gmpy.py 与 gmpy 深度集成大整数与有理数运算会直接落到 GMP 的高性能 C 实现上。二、FloatSymPy 的浮点数类型2.1 创建与自定义精度SymPy 中的浮点数都是Float类的实例定义于 sympy/core/numbers.py可以用第二个参数指定自定义精度 Float(0.1) 0.100000000000000 Float(0.1, 10) 0.1000000000 Float(0.125, 30) 0.125000000000000000000000000000 Float(0.1, 30) 0.100000000000000005551115123126最后一个例子揭示了关键事实Python 的float只有约 15 位有效数字的精度。作为输入时0.1这个 Python 浮点数的真实二进制值是0.100000000000000005551115123126...而分母为 2 的幂的数如0.125 1/8在二进制下可以精确表示因此Float(0.125, 30)显示为干净的 30 位零。相关说明在 Float 类文档字符串 中有更详细的演示如Float(0.3, 20)得到0.29999999999999998890。2.2 从字符串或Rational构造高精度Float要从十进制字符串获得真正的高精度Float应传入字符串、Rational或对Rational调用evalf Float(0.1, 30) 0.100000000000000000000000000000 Float(Rational(1, 10), 30) 0.100000000000000000000000000000 Rational(1, 10).evalf(30) 0.100000000000000000000000000000Float还支持传入空字符串自动统计有效数字位数只对字符串、int、long 有效并允许数字中出现空格或下划线 Float(123 456 789.123_456, ) 123456789.123456 Float(12e-3, ) 0.0122.3 精度如何参与运算一个数的精度决定两件事1与该数做算术运算时采用的精度2打印该数时显示的位数。当两个不同精度的数参与同一运算时结果取两者中较高的精度。官方文档给出的例子是 Float(0.1, 3)*Float(3.1415, 5) 0.31417这里0.1 /- 0.001与3.1415 /- 0.0001的乘积本身约有 0.003 的不确定度却显示了 5 位精度。因此官方明确提醒显示精度不应被当作误差传播或有效数字significance arithmetic的模型这套机制存在的目的是保证数值算法的稳定性——每次运算都以双方最高精度执行避免中间结果被截断污染。2.4 用N/evalf调整已有Float的精度 N(3.5) 3.50000000000000 N(3.5, 5) 3.5000 N(3.5, 30) 3.50000000000000000000000000000需要注意的是提升精度并不会提升准确性——底层值不变只是显示位数变多而evalf内部的确会在求值时提升工作精度但对一个本身只有 1 位精度的Float(0.1, 1).evalf(5)只能给出0.099609这种最接近真值的 5 位近似见 numbers.py 的说明。三、精度与误差处理3.1 误差跟踪与自动提精Fibonacci 案例当输入N/evalf的表达式复杂时数值误差传播成为核心问题。官方文档以第 100 个斐波那契数与近似公式 φ¹⁰⁰/√5φ 为黄金比例之差为例。用普通浮点运算相减时两数完全相同导致彻底抵消 a, b GoldenRatio**1000/sqrt(5), fibonacci(1000) float(a) 4.34665576869e208 float(b) 4.34665576869e208 float(a) - float(b) 0.0而N/evalf会跟踪误差并自动提高内部工作精度从而得到正确结果 N(fibonacci(100) - GoldenRatio**100/sqrt(5)) -5.64613129282185e-22其底层机制在 evalf_add每个操作数求值后携带各自的精度/误差信息内部表示为(re, im, re_acc, im_acc)四元组其中re_acc是相对精度的对数估计见 文件头注释complex_accuracy计算整体精度若未达到目标精度target_prec则以prec prec max(10 2**i, target_prec - acc)逐步加大工作精度重试直到达标或触及maxprec上限。3.2 精度上限与maxn数值求值**无法区分恰好为零与仅仅非常小**的表达式因此工作精度默认被限制在约 100 位十进制源码中对应 DEFAULT_MAXPREC 333即约 333 bit 二进制精度注释里还调侃道真男人把它设为 INF。尝试第 1000 个斐波那契数时 N(fibonacci(1000) - (GoldenRatio)**1000/sqrt(5)) 0.e85返回结果缺少有效数字说明N未能在默认上限内达到完全精度——只能告诉我们表达式量级小于 10⁸⁴这并不算好答案。此时可用maxn关键字强制提高工作精度上限 N(fibonacci(1000) - (GoldenRatio)**1000/sqrt(5), maxn500) -4.60123853010113e-210maxn通常可以设得很高数千位但极端情况下会显著拖慢计算。源码中maxn通过options[maxprec] max(prec, int(maxn*LG10))生效LG10 log2(10)而每次循环里的maxprec还会被进一步限制为min(oldmaxprec, 2*prec)。3.3strictTrue与PrecisionExhausted也可以设置strictTrue当达不到请求精度时抛出异常而非静默返回低精度结果 N(fibonacci(1000) - (GoldenRatio)**1000/sqrt(5), strictTrue) Traceback (most recent call last): ... PrecisionExhausted: Failed to distinguish the expression: -sqrt(5)*GoldenRatio**1000/5 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875 from zero. Try simplifying the input, using chopTrue, or providing a higher maxn for evalfPrecisionExhausted继承自ArithmeticErrorevalf.py L64由 check_target 在options[strict]为真时触发若complex_accuracy(result) prec即整体精度低于目标精度就抛出该异常并提示三种解法简化输入、使用chopTrue、提供更大的maxn。即使补全了比内公式Binet 公式的完整形式使表达式数学上精确为零N也并不知道这一点 f fibonacci(100) - (GoldenRatio**100 - (GoldenRatio-1)**100)/sqrt(5) N(f) 0.e-104 N(f, maxn1000) 0.e-13363.4chopTrue把微小量截为零在已知会发生这种抵消的场景下chop选项非常有用它把实部或虚部中非常小的数替换为精确的 0 N(f, chopTrue) 0 N(3 I*f, chopTrue) 3.00000000000000从源码看chop由 chop_parts 实现当某分量的fastloglog2 量级估计小于-prec 4时直接清零若某分量既不准确又相对很小也会被截断。chop还可以接受一个数值作为容差chopTrue时容差取标准精度否则按公式int(round(-3.321*log10(chop) 2.5))换算为精度阈值这在文档字符串中有示例 x 1e-4 N(x, chop1e-5) 0.000100000000000000 N(x, chop1e-4) 03.5 去除无意义数字round与重新求值想去除无意义数字时可以重新求值或使用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.9563.6 附evalf的完整参数表综合 EvalfMixin.evalf 签名 与官方文档evalf(n15, subsNone, maxn100, chopFalse, strictFalse, quadNone, verboseFalse)各参数含义如下参数默认值作用n15目标精度十进制位数subsNone以字典形式为符号代入数值如subs{x:3, y:1pi}maxn100允许的最大临时工作精度十进制位数chopFalse是否/以何种容差将微小实部或虚部替换为精确零strictFalse若任一子结果无法达到完全精度则抛出PrecisionExhaustedquadNone积分算法默认 tanh-sinh振荡积分用oscverboseFalse打印调试信息其中subs参数尤其值得注意官方文档在源码注释evalf L1614-L1632中特别强调用.subs()做替换可能因 Float 精度问题破坏结果例如(x y - z).subs(values)对{x:1e16, y:1, z:1e16}返回0先加 1 再减被截断而(x y - z).evalf(subsvalues)则正确返回1.00000000000000因为evalf的subs会先做evalf_subs符号替换、再在受控的高精度下求值。此外subs若以序列形式传入会抛出TypeError(subs must be given as a dictionary)。四、级数与积分的数值求值4.1 与普通闭式表达式无异求和尤其无穷级数与积分可以像普通闭式表达式一样使用并支持任意精度求值 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最后一个例子巧妙地通过高斯积分平方验证了 π 的值。从源码看Integral、Sum、Product都在evalf_table中注册了专用求值器evalf_integralL1173、evalf_sumL1331、evalf_prodL1322。它们都采用评估-检查精度-不足则提高工作精度重试的自适应循环。4.2 积分默认 tanh-sinh 算法默认情况下积分使用tanh-sinh 正交quadrature算法源码 do_integral L1141-L1143 调用ctx.quadts。该算法对光滑被积函数甚至端点奇异的积分非常高效稳健但可能难以处理高度振荡或区间内部不连续的积分。多数情况下evalf/N能正确估计误差。下面这个积分结果准确但只有 4 位有效数字 f abs(sin(x)) Integral(abs(sin(x)), (x, 0, 4)).evalf() 2.346更好的做法是把积分拆成两段被积函数在 xπ 处不连续 (Integral(f, (x, 0, pi)) Integral(f, (x, pi, 4))).evalf() 2.34635637913639类似的振荡积分例子 Integral(sin(x)/x**2, (x, 1, oo)).evalf(maxn20) 0.54.3 振荡积分quadosc上述振荡积分可以更高效地通过告诉evalf/N使用振荡正交算法解决 Integral(sin(x)/x**2, (x, 1, oo)).evalf(quadosc) 0.504067061906928 Integral(sin(x)/x**2, (x, 1, oo)).evalf(20, quadosc) 0.50406706190692837199振荡正交要求被积函数含有cos(axb)或sin(axb)因子。源码 do_integral L1127-L1140 会先用Wild模式匹配cos(A*x B)*D或sin(A*x B)*D匹配失败则抛出ValueError匹配成功后调用ctx.quadosc(f, [xlow, xhigh], periodperiod)进行带周期的振荡求积period 2*pi/A。注意许多其他形式的振荡积分可通过变量代换转化为该形式 init_printing(use_unicodeFalse) intgrl Integral(sin(1/x), (x, 0, 1)).transform(x, 1/x) intgrl oo / | | sin(x) | ------ dx | 2 | x | / 1 N(intgrl, quadosc) 0.5040670619069284.4 级数直接求和、外推与超几何快速算法无穷级数在收敛足够快时直接求和否则使用外推方法通常是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(11/k), (k, 1, oo)).evalf() 0.577215664901533 Sum(1/k-log(11/k), (k, 1, oo)).evalf(50) 0.57721566490153286060651209008240243104215933593992 EulerGamma.evalf(50) 0.57721566490153286060651209008240243104215933593992第二个例子用欧拉-马歇罗尼常数定义级数验证了EulerGamma两值精确一致到 50 位。Euler-Maclaurin 公式也用于有限级数无需逐项计算全部项即可快速逼近 Sum(1/k, (k, 10000000, 20000000)).evalf() 0.693147255559946从源码 evalf_sum L1352-L1370 可见当快速超几何路径不可用时evalf_sum会在m n 2**i * preci 从 1 到 4逐步加密的采样下调用expr.euler_maclaurin(mm, nn, epseps, eval_integralFalse)直到误差估计err eps其中eps 2**(-prec)。文档同时提醒evalf的一些默认假设并非总是最优若需要对数值求和做精细控制可直接手动调用Sum.euler_maclaurin方法。4.5 有理超几何级数的特殊优化对于有理超几何级数通项是多项式、幂、阶乘、二项式系数等的乘积存在特殊优化N/evalf能以极快速度求和到高精度。源码中对应 hypsum 函数先用hypersimp求相邻项之比并检查收敛性check_convergence若为几何级数或更快则直接整数移位求和对多项式收敛级数则用 mpmath 的nsum(..., methodrichardson)做 Richardson 外推并采用 4 倍精度工作以抵御外推过程中的抵消通过迭代加倍直到前后结果稳定。官方文档给出的例子是用拉马努金Ramanujan的 π 公式一条命令即可在不到一秒内求和到 10,000 位 f factorial n Symbol(n, integerTrue) R 9801/sqrt(8)/Sum(f(4*n)*(110326390*n)/f(n)**4/396**(4*n), ... (n, 0, oo)) N(R, 10000) #doctest: SKIP 3.141592653589793238462643383279502884197169399375105820974944592307816406286208 99862803482534211706798214808651328230664709384460955058223172535940812848111745 02841027019385211055596446229489549303819644288109756659334461284756482337867831 ...注意此处必须声明n为整数符号integerTrue超几何识别依赖这一假设。五、数值反推公式nsimplify5.1 功能与算法概览nsimplify函数定义于 sympy/simplify/simplify.py L1404尝试找出与给定输入数值相等的公式。它有两种用途对近似浮点输入猜测一个精确公式如0.1 → 1/10对复杂的符号输入猜测一个更简单的公式。其算法能够识别简单分数、简单代数表达式、给定常数的线性组合以及上述形式的某些初等函数变换。从源码看nsimplify会先将表达式在 30 位精度prec 30下求值源码注释明确要求输入应能求值到至少 30 位精度再对常数库做连分数/线性组合等搜索tolerance未指定时取输入中最不精确的值设定容差Python 浮点默认 15 位精度即tolerance10**-15。5.2 基本用法 nsimplify(0.1) 1/10 nsimplify(6.28, [pi], tolerance0.01) 2*pi nsimplify(pi, tolerance0.01) 22/7 nsimplify(pi, tolerance0.001) 355 --- 113 nsimplify(0.33333, tolerance1e-4) 1/3 nsimplify(2.0**(1/3.), tolerance0.001) 635 --- 504 nsimplify(2.0**(1/3.), tolerance0.001, fullTrue) 3 ___ \/ 2可以看出第二个参数constants允许传入可用的常数列表如[pi]tolerance设得越紧得到的公式越精确22/7→355/113fullTrue会执行更广泛的搜索从而在低容差下找到更简单的形式有理近似635/504→ 精确的2**(1/3)。5.3 高级示例官方文档给出的一组更进阶的例子 nsimplify(Float(0.130198866629986772369127970337,30), [pi, E]) 1 ---------- 5*pi ---- 2*e 7 nsimplify(cos(atan(1/3))) ____ 3*\/ 10 -------- 10 nsimplify(4/(1sqrt(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最后一组例子展示了nsimplify的反推威力给一个 30 位的Float配合常数[pi, E]能还原出1/(5*pi/7 2*E)给I**I配[pi]得到exp(-pi/2)甚至对符号级数Sum(1/n**2, (n, 1, oo))配[pi]识别出 π²/6。此外当输入含有自由符号或rationalTrue时nsimplify会退化为把浮点替换为其有理数等价形式rational_conversionbase10默认按十进制字符串表示转换exact则按精确的二进制表示转换可参看源码示例nsimplify(0.333333333333333, rationalTrue, rational_conversionexact)的结果差异。六、小结数值求值的两个入口.evalf()与N()以及别名.n()完全等价N内部先做sympify(x, rationalTrue)再调evalf默认 15 位精度n参数可指定任意位数底层由 mpmath 提供高精度算术并辅以保护位 自适应提精机制Float的精度同时决定运算精度与显示位数从字符串或Rational构造才能获得真正的高精度十进制值面对精度耗尽0.e85式结果时用maxn提高工作精度上限、strictTrue强制报错、chopTrue截断微小量或用round去尾积分默认走 tanh-sinh 正交振荡/间断被积函数可拆分区间或改用quadosc级数求值综合运用直接求和、Euler-Maclaurin、Richardson 外推与超几何快速算法配合Sum.euler_maclaurin可做精细控制nsimplify是数值反推符号公式的利器配合constants、tolerance、full参数可在工程与教学场景中快速验证猜想。更深入的实现细节可直接查阅 sympy/core/evalf.py自适应求值核心、evalf_table分派表、hypsum/do_integral/evalf_sum等与 sympy/core/numbers.pyFloat类并在 sympy/core/tests/test_evalf.py 等测试文件中找到大量可复现的求值用例。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考