数值分析的Sage快速入门¶

这个 Sage 快速入门教程是为MAA准备工作坊“Sage：对本科生使用开源数学软件”开发的（由NSF提供资金，到期日：0817071）。它是根据Creative Commons Attribution -ShareLiked 3.0许可证授权的 (CC BY-SA ）

Sage包括许多用于数值分析研究的工具。

开始位置是Sage中默认的十进制数类型。

基本分析¶

这个 RealField 使用任意精度初始化（用 MPFR ）
默认实数 (RR 是 RealField(53) （即53位精度）。
但你可以让他们活在你想要的任何精确的地方。

sage: ring=RealField(3)

要打印实际数字（不舍入最后几个不精确的数字只显示正确的数字），请调用 .str() 方法：

sage: print(ring('1').nextabove())
1.2

sage: print(ring('1').nextabove().str())
1.2
sage: print(ring('1').nextbelow().str())
0.88

让我们改变我们的精度。

sage: ring=RealField(20)
sage: print(ring('1').nextabove().str())
1.0000019
sage: print(ring('1').nextbelow().str())
0.99999905

也可以指定舍入模式。

sage: ringup=RealField(3,rnd='RNDU')
sage: ringdown=RealField(3,rnd='RNDD')

sage: ring(1/9).str()
'0.11111116'

sage: ringup(1/9).str()
'0.13'

sage: ringdown(1/9).str()
'0.10'

sage: ring(1/9).str(base=2)
'0.00011100011100011100100'

让我们看看分数到底是多少。

sage: ring(1/9).exact_rational()
233017/2097152

sage: n(233017/2097152)
0.111111164093018

转换为浮点二进制（IEEE格式）¶

sage: x=13.28125

给你 x 作为二进制浮点数。

sage: x.str(base=2)
'1101.0100100000000000000000000000000000000000000000000'

从这个二进制浮点表示，我们可以构造 x 再一次。

sage: xbin=2^3 +2^2 + 2^0+2^-2+2^-5; xbin
425/32

sage: xbin.n()
13.2812500000000

为了检查，让我们问塞奇什么 x 是一个精确的分数（不涉及四舍五入）。

sage: x.exact_rational()
425/32

现在让我们转换 x 计算机中常用的ieee754二进制格式。为 IEEE 754 ，获得二进制格式的第一步是规范化数字，或者将数字表示为1到2之间的数字乘以2的幂。为了我们 x 上面，我们乘以 2^{{-3}} 得到一个介于1和2之间的数字。

sage: (x*2^(-3)).str(base=2)
'1.1010100100000000000000000000000000000000000000000000'

标准化数的小数部分称为尾数（或） 有效数 ）

sage: f=(x*2^(-3)).frac() # .frac() gets the fraction part, the part after the decimal
sage: mantissa=f.str(base=2)
sage: mantissa
'0.10101001000000000000000000000000000000000000000000000'

因为我们乘以 \(2^{{-3}}\) 为了规范化数字，我们需要存储这些信息。为了不必担心存储负数，我们添加了1023。在下面， c 将成为我们的代表。

sage: c=(3+1023).str(base=2)
sage: c
'10000000010'

注意 c 有11位，这正是我们想要的。

sage: len(c)
11

评价 mantissa[2:54] 将给出尾数小数点后的前52位二进制数字。注意，我们不需要在小数点之前存储前导1，因为从我们规范化事物的方式来看，它总是在那里。这使得我们仅使用52位存储器就可以获得53位精度。

sage: len(mantissa[2:54])
52

因为原来的数字是正数，我们的符号位是零。

sage: sign='0'

这里是我们的64位双精度浮点数。

sage: sign+' '+c+' '+mantissa[2:54] # the [2:] just chops off the '0.', since we just need to store the digits after the decimal point
'0 10000000010 1010100100000000000000000000000000000000000000000000'

sage: len(sign+c+mantissa[2:54]) # it's 64 bits!
64

在这里，我们将从我们构造的浮点表示转换回原始数字。

sage: ((-1)^(int(sign)) * 2^(int(c,base=2)-1023)*(1+RR(mantissa[:54], base=2)))
13.2812500000000

sage: x
13.2812500000000

所以他们同意了！

Sage使用先进的数字库MPFR，使用用户指定的任何精度执行精确的浮点运算。MPFR对于规范化有一个稍微不同的约定。在MPFR中，我们通过将尾数乘以适当的2的幂，使尾数成为整数，而不是二进制分数。这允许我们使用大整数库和复杂的技术以任意精度进行计算。

sage: x.sign_mantissa_exponent()
(1, 7476679068876800, -49)

sage: 7476679068876800*2^(-49)
425/32

注意这里的尾数和上面的尾数有相同的零位/非零位（在我们切掉前面的1之前）。

sage: 7476679068876800.str(base=2)
'11010100100000000000000000000000000000000000000000000'

区间算术¶

Sage还允许您使用间隔来计算，以跟踪错误边界。这些基本上使用上面显示的向上取整和向下取整功能。

sage: ring=RealIntervalField(10)
sage: a=ring(1/9)
sage: a
0.112?

问号表示法是指数字包含在通过递增和递减数字的最后一位而得到的间隔中。见 documentation for real interval fields 有关详细信息。在上面的例子中，Sage说1/9在0.111和0.113之间。下面，我们看到了 1/a 在8.9到9.1之间。

sage: 1/a
9.0?

如果我们显式地打印出区间，我们可以得到更精确的区间估计。

sage: print((1/a).str(style='brackets'))
[8.9843 .. 9.0157]

包含的软件¶

Scipy（包含在Sage中）有很多数值算法。看到了吗 the Scipy docs .

Mpmath也包含在Sage中，并且包含了大量的数字内容。看到了吗 the mpmath codebase .

这个 Decimal python module 对于课本练习也很有用，这些练习涉及到以10为基数的四舍五入。

精确绘图¶

有时，绘图涉及一些相当糟糕的舍入错误，因为绘图计算是用机器精度的浮点数完成的。

sage: f(x)=x^2*(sqrt(x^4+16)-x^2)
sage: plot(f,(x,0,2e4))
Graphics object consisting of 1 graphics primitive

我们可以使用 fast_callable 系统。

sage: R=RealField(100) # 100 bits
sage: g=fast_callable(f, vars=[x], domain=R)
sage: plot(g,(x,0,2e4))
Graphics object consisting of 1 graphics primitive