【数值计算库-超长笔记】Python-Mpmath库：高精度数值计算

原文链接：https://www.cnblogs.com/aksoam/p/18279394
更多精彩，关注博客园主页，不断学习！不断进步！

我的主页

csdn很少看私信，有事请b站私信

博客园主页-发文字笔记-常用

有限元鹰的主页 内容：

1. ABAQUS数值模拟相关
2. Python科学计算
3. 开源框架，编程学习笔记

哔哩哔哩主页-发视频-常用

FE-有限元鹰的个人空间 内容：

1. 模拟案例
2. 网格划分
3. 游戏视频，及其他搬运视频

Numpy 和 MpMath 都是 Python 中常用的数学库，它们都提供了许多数学函数和工具，但它们的设计目的和实现方式有所不同。Numpy 主要用于处理大规模的数值数据，而 MpMath 则专注于高精度计算和符号计算。

虽然 Numpy 和 MpMath 的设计目的不同，但它们之间可以进行互操作。具体来说，可以将 Numpy 数组转换为 MpMath 的高精度数值类型，然后使用 MpMath 的函数进行计算。同样地，也可以将 MpMath 的数值类型转换为 Numpy 数组，然后使用 Numpy 的函数进行计算。

Python-Mpmath库：高精度数值计算

Checking that it works

# import mpmath as mp
from mpmath import *
from icecream import ic
mp.dps = 50         # Set precision
print(mpf(2)**mpf('0.5'))

1. Introduction

1.2 Basic usage of mpmath

mpmath的数值类型: mpf(real float),mpc(complex float),matrix

print(mp.dps)
print(mp.cos(3.14))a1=mpf(4)
a2=mpf('4.0')
ic(a1,a2)
a3=mpf("1.25e6")
ic(a3)
a4=mpf("-inf")
ic(a4)a5=mpc(1,2)
ic(a5)
a6=mpc(complex(1,20))
ic(a6)ic(a5.conjugate())  # 共轭复数
ic(a5.real)         # 实部
ic(a5.imag)         # 虚部

Mpmath使用全局工作精度;它不跟踪单个数字的精度或准确性。执行算术运算或调用mpf()将结果舍入到当前的工作精度。工作精度由一个名为mp的上下文对象控制

print(mp)
# Mpmath settings:
#   mp.prec = 169               [default: 53]
#   mp.dps = 50                 [default: 15]
#   mp.trap_complex = False     [default: False]# prec---二进制精度
# dps---十进制精度
#  prec = 3.33*dps# 当创建一个新的mpf时，该值最多将与输入一样准确
#!!警告： 将mpmath数字与Python浮点数混合使用时要小心
mp.dps = 30
# bad
ic(mpf(10.8))# good
ic(mpf('10.8'))# also good
ic(mpf(108)/mpf(10))# ? 然而，0.5、1.5、0.75、0.125等通常作为输入是安全的，因为它们可以用Python浮点数精确地表示

1.3 输出格式化

# 设置 mp.pretty = True 可以输出更好看的结果mp.pretty = True
mpf(0.6)
# output: 0.599999999999999977795539507497mp.pretty = False
mpf(0.6)
# output: mpf('0.599999999999999977795539507496869')

1.4 输出的小数点位数

use mpmath.nstr() and mpmath.nprint() 函数.默认情况下打印数字的位数由工作精度决定。要指定要显示的位数，而不改变工作精度，

a=mpf(1)/6.0mp.nstr(a, 5)
# output: '0.16667'mp.nprint(a, 10)
# output: '0.1666666667'

2. BASIC FEATURES

mpmath中的高级代码是作为“上下文对象”上的方法实现的。上下文实现算术、类型转换和其他基本操作。上下文还保存诸如精度之类的设置，并存储缓存数据。提供了一些不同的上下文(具有基本兼容的接口)，以便高级算法可以与底层算法的不同实现一起使用，从而允许不同的特性和速度-精度权衡。

目前，mpmath提供了以下上下文:

• Arbitrary-precision arithmetic (mp)
• A faster Cython-based version of mp (used by default in Sage, and currently only available there)
• Arbitrary-precision interval arithmetic (iv)
• Double-precision arithmetic using Python’s builtin float and complex types (fp)

2.1 Arbitrary-precision floating-point (mp)

mp支持任意精度,并且支持绝大多数函数,经过充分测试,充分优化.

2.2 Arbitrary-precision interval arithmetic (iv)

iv支持区间运算,但是在版本1.3.0中,iv还处于测试阶段,只支持少部分函数.

2.3 Double-precision arithmetic using Python’s builtin float and complex types (fp)

尽管mpmath通常是为任意精度算术设计的，但许多高级算法可以很好地处理普通的Python浮点数和复数，它们使用硬件双精度(在大多数系统上，这相当于53位精度).mp object的函数会把输入转换为mpmath numbers type. fp object则是将输入转换为Python的浮点数和复数类型.当需要进行大量的函数求值(数值积分、绘图等)，并且当fp算法提供足够的精度时，这可以大大提高mp算法的速度。

基于以上优势,可以使用fp.func 进行运算.(使用func 或mp.func 都是mp的函数)

使用fp算法计算的结果可能比使用等效精度(mp)的mp计算的结果不那么精确.(mp.prec = 53)，因为后者经常使用更高的内部精度。其精度高度依赖于:

• for some functions, fp almost always gives 14-15 correct digits; for others, results can be accurate to only 2-3 digits or even completely wrong.
• The recommended use for fp is therefore to speed up large-scale computations where accuracy can be verified in advance on a subset of the input set, or where results can be verified afterwards.

fp.pretty=True
u=fp.rand()print(u)
# 0.6473090346719316print(type(u))
# <class 'float'>m=fp.matrix([[1,2],[3,4]])**2print(m)
# [ 7.0  10.0]
# [15.0  22.0]print(type(m))
# <class 'mpmath.matrices.matrices.matrix'>

2.4 基本函数

fp bject 为基本函数 包装了Python的math和cmath模块。支持实数和复数，并自动为实数输入生成复数结果(数学会引发异常)

fp.sqrt(5)
# 2.23606797749979fp.sqrt(-5)
# 2.23606797749979jfp.sin(fp.pi/2)
# 1.0fp.power(-2,0.5)
# (8.659560562354934e-17+1.4142135623730951j)print(type(fp.power(-2,0.5)))
# <class 'complex'>

2.4 实用函数

2.4.1 converion and printing

# mp.mpmathify(x,strings=True)
# 能够将x转换为mpf or mpc类型fp.dps=15
fp.pretty = Truempmathify(3.5)
# mpf('3.5')mpmathify(3+4j)
#mpc(real='3.0', imag='4.0')mpmathify('3.5')
# mpf('3.5')# mpmath.nstr(x,n=6,**kwargs)
# 将mpf或mpc转换为具有n位有效数字的十进制字符串
# nprint(x,n=6) 会直接打印nstr()的字符串结果f=mp.mpf(1.53)
print(f)
nprint(f,n=12)
# 1.530000000000000026645352591
# 1.53

2.4.2 算术运算

包括加减乘除、取模、幂运算等。

# ! mpmath.fadd(x, y,**kwargs)
# 实现x+y,返回浮点数结果.可以自定义精度和舍入模式.默认精度是工作精度,
# kwargs中可以设置精度和舍入模式.:
# + prec,dps(int)--设置精度; exact=True--不四舍五入.
# + rounding(str)--设置舍入模式,可选'n'(default)--nearst,'f'--floor,'c'-ceiling,'d'-down,'u'-up.
x=mp.rand()
y=mp.rand()
ic(x,y)
nprint(mp.fadd(x,y,dps=5,rounding='u'))
# ic| x: mpf('0.936352251244129142275205120579475')
#     y: mpf('0.118268786487759443295747691143063')
# 1.05462# ! mpmath.fsub(x, y,**kwargs)
# 实现x-y,返回浮点数结果.
x=mp.mpf(1.53)
y=mp.mpf(2.87)
print(x-y)
# -1.34000000000000007993605777301
print(mp.fsub(x,y))
# -1.34000000000000007993605777301# ! mpmath.fneg(x,**kwargs)
# 实现-x,返回浮点数结果.
print(mp.fneg(mp.mpf(2.5)))
# -2.5# ! mpmath.fmul(x, y,**kwargs)
# 实现x*y,返回浮点数结果.
x=mp.mpf(2.5)
y=mp.mpf(3.14)
print(x*y)
# 7.85000000000000031086244689504
print(mp.fmul(x,y))
# 7.85000000000000031086244689504a=mp.mpc(2,5)
b=mp.mpc(3,4)
print(a*b)
# (-14.0 + 23.0j)
print(mp.fmul(a,b))
# (-14.0 + 23.0j)# ! mpmath.fdiv(x, y,**kwargs)
# 实现x/y,返回浮点数结果.
x=mp.mpf(3)
y=mp.mpf(2)
print(x/y)
# 1.5
print(mp.fdiv(x,y))
# 1.5# ! mpmath.fmod(x,y)
# Converts 𝑦 and 𝑧 to mpmath numbers and returns 𝑦 mod 𝑧. 
# For mpmath numbers, this is equivalent to x % y.
mp.fmod(100,3)
# mpf('1.0')
print(100%3)
# 1# ! mpmath.fsum(terms,absolute=False,squared=False,**kwargs)
# 计算一个array-like的terms的和.如果len(terms)>2,这会比python内置的sum函数快很多.
# + absolute=True--返回每一项的绝对值的总和.
# + squared=True--返回每一项的平方的总和.
mp.fsum([1,2,3,4,5])
# mpf('15.0')mp.fsum([1,-2,3,-4,5],absolute=True)
# mpf('15.0')
mp.fsum([1,-2,3,-4,5],absolute=False)
# mpf('3.0')mp.fsum([1,2,3,4,5],squared=True)
# mpf('55.0')# ! mpmath.fprod(factors)
# 计算一个array-like object的每一项的连续乘积.
# 如: fprod([1,2,3,4,5])=1*2*3*4*5=120
mp.fprod([1,2,3,4,5])
# mpf('120.0')# ! mpmath.fdot(A, B=None, conjugate=False)
# 接受两个array-like object A和B,计算A,B的点积. sum(A[i]*B[i]) ,i=0,...,len(A)-1.
# if conjugate is True, B is conjugated before multiplication.  
# 如: fdot([1,2,3],[4,5,6])=1*4+2*5+3*6=32
mp.fdot([1,2,3],[4,5,6])
# mpf('32.0')

2.4.3 复数相关函数

fabs(x) :返回复数x的模或者实数x的绝对值.
arg(x):返回x的相位角(以弧度表示).-pi< arg(x) < pi .
sign(x) :返回x的符号,sign(x)=x/|x|.
re(x) :返回x的实部.
im(x) :返回x的虚部.
conj(x) :返回x的共轭复数$\bar{x}$
polar(x) :将复数x分解为极坐标形式(r,θ)$x=r\exp (i\theta )$
rect(r,θ) :将极坐标形式(r,θ)转换为复数形式$x=r(\cos \theta +i\sin \theta )$

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> chop(rect(2, pi))
-2.0
>>> rect(sqrt(2), -pi/4)
(1.0 - 1.0j)

2.4.4 处理整数和小数

floor(x) : 向下取整数
ceil(x) : 向上取整数
nint(x) : 四舍五入取整数
frac(x) : 返回x的小数部分,frac(x)=x-floor(x)

>>> from mpmath import *
>>> mp.pretty = False
>>> frac(1.25)
mpf('0.25')
>>> frac(3)
mpf('0.0')
>>> frac(-1.25)
mpf('0.75')

2.4.5 公差和近似比较

chop(x, tol=None) : 去除小的实部或者虚部,或者将接近于0的转换为精确的0,x可以是数字或者iterable.The tolerance defaults to 100*eps.
almosteq(s, t, rel_eps=None, abs_eps=None) : 判断s-t是否小于给定的误差,rel_eps和abs_eps分别是相对误差和绝对误差.绝对误差=|s-t|,相对误差=|s-t|/max(|s|,|t|).

2.4.6 数组生成

fraction(p,q) ;给定两个整数p和q,返回分数p/q.比其他的精度更好

>>> from mpmath import *
>>> mp.dps = 15
>>> a = fraction(1,100)
>>> b = mpf(1)/100
>>> print(a); print(b)
0.01
0.01
>>> mp.dps = 30
>>> print(a); print(b) # a will be accurate
0.01
0.0100000000000000002081668171172
>>> mp.dps = 15

rand() : 随机生成一个mpf随机数,范围在[0,1)之间.
arrange(*args):返回mpf数组序列. 最后一个数不是b

• arrange(b):[0,1,2,…,x]
• arrange(a,b):[a,a+1,a+2,…,x]

$b-1\le x<b$

• arrange(a,b,h):[a,a+h,a+2h,…,x]

$b-h\le x<b$

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> arange(4)
[mpf('0.0'), mpf('1.0'), mpf('2.0'), mpf('3.0')]
>>> arange(1, 2, 0.25)
[mpf('1.0'), mpf('1.25'), mpf('1.5'), mpf('1.75')]
>>> arange(1, -1, -0.75)
[mpf('1.0'), mpf('0.25'), mpf('-0.5')]

linspace(start,end,num) :返回浮点数数组,从start到end均匀分成num个元素.默认包括end.元素之间的间隔$\Delta =(end-start)/(num-1)$

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> linspace(1, 4, 4)
[mpf('1.0'), mpf('2.0'), mpf('3.0'), mpf('4.0')]

使用 endpoint=False 可以排除end.此时: 元素之间的间隔$\Delta =(end-start)/(num+1)$

>>> linspace(1, 4, 4, endpoint=False)
[mpf('1.0'), mpf('1.75'), mpf('2.5'), mpf('3.25')]

2.5 plotting函数

plot([cos, sin], [-4, 4]): 绘制cos和sin函数的图像,x轴范围为[-4,4].

2.5.1 2D plotting

This function requires matplotlib (pylab).

plot(f, xlim=[- 5, 5], ylim=None, points=200, file=None, dpi=None, singularities=[], axes=None)

给定一个函数f(x),或者一个函数序列[f1(x),f2(x),…,fn(x)],并指定x周范围xlim=[start,end],绘制2D图像

使用singularities可以在指定点,消除数值奇异或者无穷大值.例如tan函数

# 使用singularities绘制tan函数
from mpmath import *
plot(tan, [-2*pi, 2*pi],singularities=[-pi/2, pi/2,-1.5*pi,1.5*pi],points=10000)

# 不使用singularities绘制tan函数
from mpmath import *
plot(tan, [-2*pi, 2*pi],points=10000)

例如:

plot(lambda x: exp(x)*li(x), [1, 4])
plot([fresnels, fresnelc], [-4, 4])
plot([sqrt, cbrt], [-4, 4])
plot(lambda t: zeta(0.5+t*j), [-20, 20])
plot([floor, ceil, abs, sign], [-5, 5])

3. 高级特性

3.1 矩阵

3.1.1 创建矩阵

mpmath中的矩阵是使用字典实现的。只存储非零值，因此表示稀疏矩阵的成本很低。

最简单的矩阵创建方法是使用matrix()函数,指定矩阵维度,可以得到空矩阵;也可以基于列表创建矩阵,这在根据已知numpy数组创建矩阵时很有用.

# 创建空矩阵
>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> matrix(2)
matrix(
[['0.0', '0.0'],
['0.0', '0.0']])
>>> matrix(2, 3)
matrix(
[['0.0', '0.0', '0.0'],
['0.0', '0.0', '0.0']])
# 基于列表创建矩阵
>>> matrix([[1, 2], [3, 4]])
matrix(
[['1.0', '2.0'],
['3.0', '4.0']])
>>> mp.matrix(np.array([[1,2,5],[3,4,9]]))
matrix(
[['1.0', '2.0', '5.0'],['3.0', '4.0', '9.0']])
>>>mp.matrix(np.array([1,2,5,6]))
matrix(
[['1.0'],['2.0'],['5.0'],['6.0']])
>>>mp.matrix(np.array([[1,2,5,6]]))
matrix(
[['1.0', '2.0', '5.0', '6.0']])

访问matrix object的rows,cols属性可以得到矩阵的行数,列数.直接修改rows,cols可以改变矩阵的维度.(新增的元素为0,删除的元素会丢失)

>>> A = matrix(3, 2)
>>> A
matrix(
[['0.0', '0.0'],
['0.0', '0.0'],
['0.0', '0.0']])
>>> A.rows
3
>>> A.cols
2
>>> A.rows = 2
>>> A
matrix(
[['0.0', '0.0'],
['0.0', '0.0']])

索引矩阵元素的语法为:matrix[i,j],i表示行,j表示列.

>>> A = matrix(2)
>>> A[1,1] = 1 + 1j
>>> print(A)
[0.0 0.0]
[0.0 (1.0 + 1.0j)]

• 高级方法

zeros(rows, cols) : 创建rows行cols列的零矩阵
ones(rows, cols) : 创建rows行cols列的单位矩阵
eye(rows, cols) : 创建rows行cols列的单位矩阵
randmatrix(rows, cols) : 创建rows行cols列的随机矩阵

3.1.2 向量

向量也可以使用matrix表示(rows=1 or cols=1),行，列向量的创建方法:

# 创建列向量
>>> from mpmath import *
>>> a=matrix([1,2,3])
>>> a
matrix(
[['1.0'],['2.0'],['3.0']])
>>> a.rows
3
>>> a.cols
1
# 创建行向量
>>> b=matrix([[1,2,3]])
>>> b
matrix(
[['1.0', '2.0', '3.0']])
>>> b.rows
1
>>> b.cols
3

3.1.3 矩阵复制

• 直接打印矩阵,控制有效位数:

>>> print(randmatrix(3))
[ 0.782963853573023 0.802057689719883 0.427895717335467]
[0.0541876859348597 0.708243266653103 0.615134039977379]
[ 0.856151514955773 0.544759264818486 0.686210904770947]
>>> nprint(randmatrix(5), 3)
[2.07e-1 1.66e-1 5.06e-1 1.89e-1 8.29e-1]
[6.62e-1 6.55e-1 4.47e-1 4.82e-1 2.06e-2]
[4.33e-1 7.75e-1 6.93e-2 2.86e-1 5.71e-1]
[1.01e-1 2.53e-1 6.13e-1 3.32e-1 2.59e-1]
[1.56e-1 7.27e-2 6.05e-1 6.67e-2 2.79e-1]

由于matrix object属于mutable,因此可以直接修改元素的值,或者使用copy()方法复制矩阵.涉及到浅拷贝,深拷贝,赋值的概念,可以参考python的相关知识.

# 拷贝矩阵
>>> a=matrix([[1,2],[3,4]])
>>> a
matrix(
[['1.0', '2.0'],['3.0', '4.0']])
>>> b=a.copy()
>>> b
matrix(
[['1.0', '2.0'],['3.0', '4.0']])
>>> a[0,0]=0
>>> a
matrix(
[['0.0', '2.0'],['3.0', '4.0']])
>>> b
matrix(
[['1.0', '2.0'],['3.0', '4.0']])
# 直接赋值
>>> id(a)
2441463350224
>>> b=a
>>> id(b)
2441463350224   # 两个变量指向同一块内存空间
>>> b
matrix(
[['1.0', '2.0'],['3.0', '4.0']])
>>> b[1,1]=100
>>> a
matrix(
[['1.0', '2.0'],['3.0', '100.0']])
>>> b
matrix(
[['1.0', '2.0'],['3.0', '100.0']])

matrix object可以使用.tolist()方法转换为嵌套列表

>>> a
matrix(
[['0.0', '2.0'],['3.0', '100.0']])
>>> a.tolist()
[[mpf('0.0'), mpf('2.0')], [mpf('3.0'), mpf('100.0')]]
>>> import numpy as np
>>> np.array(a.tolist())
array([[mpf('0.0'), mpf('2.0')],[mpf('3.0'), mpf('100.0')]], dtype=object)
>>> np.array(a.tolist())[1,1]
mpf('100.0')

3.1.4 矩阵运算

• 矩阵和实数之间的加减乘除,次幂运算(A^n)

>>> A = matrix([[1, 2], [3, 4]])
>>> B = matrix([[-2, 4], [5, 9]])
# 矩阵相加
>>> A + B
matrix(
[['-1.0', '6.0'],
['8.0', '13.0']])
>>> A - B
matrix(
[['3.0', '-2.0'],
['-2.0', '-5.0']])
# 举着相加需要注意维度的一致性
>>> A + ones(3)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "...", line 238, in __add__
raise ValueError('incompatible dimensions for addition')
ValueError: incompatible dimensions for addition
# 矩阵乘法
>>> A * 2
matrix(
[['2.0', '4.0'],
['6.0', '8.0']])
>>> A / 4
matrix(
[['0.25', '0.5'],
['0.75', '1.0']])
>>> A - 1
matrix(
[['0.0', '1.0'],
['2.0', '3.0']])

• 矩阵之间的加减乘除运算,维度必须一致,否则会报错.

>>> A * B
matrix(
[['8.0', '22.0'],
['14.0', '48.0']])
# 行向量*列向量
>>> matrix([[1, 2, 3]]) * matrix([[-6], [7], [-2]])
matrix(
[['2.0']])
>>> A**2
matrix(
[['7.0', '10.0'],
['15.0', '22.0']])

• 矩阵求逆:A**(-1)

>>> A
matrix(
[['1.0', '2.0'],['3.0', '4.0']])
>>> A**-1
matrix(
[['-2.0', '1.0'],['1.5', '-0.5']])
>>> A*A**-1
matrix(
[['1.0', '1.0842021724855e-19'],['-2.16840434497101e-19', '1.0']])

• 矩阵转置:A.T

>>> A
matrix(
[['1.0', '2.0'],['3.0', '4.0']])
>>> A.T
matrix(
[['1.0', '3.0'],['2.0', '4.0']])

• 矩阵行列式

对于一个方阵,mpmath.det(A)可以计算其行列式.

>>> from mpmath import mp
>>> A = mp.matrix([[7, 2], [1.5, 3]])
>>> print(mp.det(A))
18.0

• 计算范数

矩阵和向量无法直接比较大小,但可以用范数(norm)来衡量矩阵和向量的大小.

范数:将矩阵或向量映射到非负实数值,使得该值表示矩阵或向量的大小.

mpmath.norm(x,p=2):给定x(iterable),计算x的$p范数=({\sum }_k|x_k{|}^p{)}^{1/p},1\le p\le \infty$.

特别的:如果x不是iterable,则返回absmax(x).

p=1: 绝对值范数,即${\sum }_k|x_k|$

p=2: 二范数(标准欧几里得向量范数),即$\sqrt{{\sum }_k|x_k{|}^2}$

p=inf: 最大范数,即$\max |x_k|$

需要注意的是: 如果x是矩阵,且p=2,则返回Frobenius范数,即$\sqrt{{\sum }_{i,j}|A_{ij}{|}^2}$. 对于算子矩阵规范，使用mnorm()代替

>>> A
matrix(
[['1.0', '2.0'],['3.0', '4.0']])
>>> 1+4+9+16
30
>>> sqrt(30)
mpf('5.4772255750516612')
>>> norm(A,p=2)
mpf('5.4772255750516612')

mpmath.mnorm(A, p=1): 计算矩阵A的矩阵(算子)p-范数.1.3.0版本支持p=1,inf.

特别的:

p=1: 矩阵(算子)max column sum 范数,即$max{\sum }_i|A_{ij}|$.

p=inf: 矩阵(算子)max row sum 范数,即$max{\sum }_j|A_{ij}|$.

p=‘f’: 矩阵(算子)Frobenius范数,相当于elementwise 2-norm.Frobenius范数是谱范数的近似值,且满足:

Frobenius范数缺乏一些范数应有的数学性质

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> A = matrix([[1, -1000], [100, 50]])
>>> mnorm(A, 1)
mpf('1050.0')
>>> mnorm(A, inf)
mpf('1001.0')
>>> mnorm(A, 'F')
mpf('1006.2310867787777')
# 矩阵(算子)Frobenius范数与elementwise 2-norm的比较
>>> A
matrix(
[['1.0', '2.0'],['3.0', '4.0']])
>>> norm(A,p=2)
mpf('5.4772255750516612')
>>> mnorm(A,p='f')
mpf('5.4772255750516612')

未完待续…2024.07.02