用于数值分析的SAGE快速入门¶

这 Sage 快速入门教程是为MAA预备研讨会“SAGE：在本科生中使用开源数学软件”开发的(由美国国家科学基金会提供，截止日期0817071)。它是在知识共享署名-Share Alike 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

为了检查，让我们问Sage什么 x 是一个精确的分数(不涉及舍入)。

sage: x.exact_rational()
425/32

现在让我们转换一下 x 转换为计算机中常用的IEEE 754二进制格式。为 IEEE 754 ，获取二进制格式的第一步是将数字标准化，或将数字表示为1和2之间的数字乘以2的幂。 x 在上面，我们乘以 2^{-3} 若要获取介于1和2之间的数字，请执行以下操作。

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

归一化数的分数部分称为 mantissa (或 significand )。

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 。

Mpath也包含在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

方法专门计算所有中间步骤的精度到100位。 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