mpmath

Multiple Precision math

Пакет для работы с числами с плавающей точкой со сколь угодно высокой точностью. В нём реализованы алгоритмы вычисления элементарных функций, а также большого количества специальных функций.

In [1]:
from mpmath import *
%matplotlib inline

Точность контролируется глобальным объектом mp.

In [2]:
print(mp)

Mpmath settings:
  mp.prec = 53                [default: 53]
  mp.dps = 15                 [default: 15]
  mp.trap_complex = False     [default: False]

prec - число бит в мантиссе, dps - число значащих десятичных цифр. Если изменить один из этих атрибутов, другой изменится соответственно.

In [3]:
mp.dps=50
print(mp)

Mpmath settings:
  mp.prec = 169               [default: 53]
  mp.dps = 50                 [default: 15]
  mp.trap_complex = False     [default: False]

mpf создаёт число с плавающей (multiple precision float) точкой из строки или числа.

In [4]:
x=mpf('0.1')
print(x)

0.1
In [5]:
x
Out[5]:
mpf('0.10000000000000000000000000000000000000000000000000007')

А вот так делать не надо. 0.1 сначала преобразуется в число с плавающей точкой со стандартной (т.е. двойной) точностью, а потом уже оно преобразуется в mpf.

In [6]:
y=mpf(0.1)
print(y)

0.1000000000000000055511151231257827021181583404541

Чтобы не потерять точность, нужно делать mpf из строки или из отношения целых чисел (вероятно, со знаменателем вида $10^n$).

In [7]:
y=mpf(1)/10
print(y)

0.1

Математические константы типа $\pi$ или $e$ реализованы в виде ленивых объектов. Они содержат сколько-то вычисленных бит плюс алгоритм, позволяющий получить больше бит, если потребуется.

In [8]:
pi
Out[8]:
<pi: 3.14159~>
In [9]:
mp.prec=53
print(pi)

3.14159265358979
In [10]:
mp.prec=169
print(pi)

3.1415926535897932384626433832795028841971693993751

Когда объект pi встречается в выражении, из него делается число с текущей точностью.

In [11]:
+pi
Out[11]:
mpf('3.1415926535897932384626433832795028841971693993751068')

Реализованы арифметические операции и элементарные функции.

In [12]:
sin(pi/4)**2
Out[12]:
mpf('0.50000000000000000000000000000000000000000000000000134')

Специальные функции

In [13]:
plot([lambda x:besselj(0,x),
      lambda x:besselj(1,x),
      lambda x:besselj(2,x)],[0,10])
In [14]:
plot([lambda x:legendre(0,x),
      lambda x:legendre(1,x),
      lambda x:legendre(2,x),
      lambda x:legendre(3,x)],[-1,1])
In [15]:
plot([lambda x:polylog(2,x),
      lambda x:polylog(3,x),
      lambda x:polylog(4,x)],[-4,1])
In [16]:
print(pi**2/zeta(2),pi**4/zeta(4))

6.0 90.0
In [17]:
splot(lambda x,y:abs(gamma(x+j*y)),[-4,4],[-4,4])
In [18]:
gamma(1.5)/sqrt(pi)
Out[18]:
mpf('0.5')

Решение уравнений

Корни многочлена

In [19]:
l=[1,0,0,0,1,1]
r=polyroots(l)
for x in r:
    print(x)

-0.75487766624669276004950889635852869189460661777279
(0.8774388331233463800247544481792643459473033088864 - 0.74486176661974423659317042860439236724016308490682j)
(0.8774388331233463800247544481792643459473033088864 + 0.74486176661974423659317042860439236724016308490682j)
(-0.5 + 0.86602540378443864676372317075293618347140262690519j)
(-0.5 - 0.86602540378443864676372317075293618347140262690519j)
In [20]:
for x in r:
    print(polyval(l,x))

0.0
(2.672764710092195646140536467151481878815196880105e-51 - 4.6512209026900071036135543450317217153701063147101e-51j)
(2.672764710092195646140536467151481878815196880105e-51 + 4.6512209026900071036135543450317217153701063147101e-51j)
(0.0 + 0.0j)
(0.0 + 0.0j)

Решение уравнения

In [21]:
def f(x):
    return exp(-x)-sin(x)
In [22]:
plot(f,[0,pi])
In [23]:
findroot(f,(0.5,0.7))
Out[23]:
mpf('0.58853274398186107743245204570290368853127151610903053')

Решение системы уравнений

In [24]:
findroot([lambda x,y:x**2+y**2-1,lambda x,y:x*y-1/4],(1,0.25))
Out[24]:
matrix(
[['0.9659258262890682867497431997288973676339048390084'],
 ['0.25881904510252076234889883762404832834906890131993']])

Производные

In [25]:
diff(f,0.5)
Out[25]:
mpf('-1.4841132216030061397200811175950101054335633325969314')
In [26]:
diff(f,0.5,2)
Out[26]:
mpf('1.085956198316836423877087470206751841523721503427788')
In [27]:
diff(lambda x,y:sin(x)*cos(y),(pi,pi),(1,2))
Out[27]:
mpf('-1.0')

Интегралы

При вычислении этого интеграла все вычисления будут производиться с точностью, на 5 значащих цифр большей; затем она вернётся к прежней.

In [28]:
with extradps(5):
    I=quad(lambda x:log(x)**2/(1+x),(0,1))
    print(I)

1.803085354739391428099607242267174986147479438510748322

Допустим, у нас есть причины подозревать, что этот интеграл равен $\zeta(3)$, умноженному на рациональное число (с не очень большими числителем и знаменателем). pslq([x1,x2,...]) находит целые числа $n_1$, $n_2$, ... такие, что $n_1\,x_1 + n_2\,x_2 + \cdots = 0$. Это - метод нахождения тождеств, называемый экспериментальной математикой. Для этого часто требуются вычисления с очень высокой точностью.

In [29]:
pslq([I,zeta(3)])
Out[29]:
[-2, 3]

То есть наш интеграл равен $\frac{3}{2} \zeta(3)$. Это, конечно, не доказательство. Но если мы ещё увеличим точность вычисления интеграла, а результат pslq не изменится, то мы можем быть практически уверены, что этот результат верен.

Двойной интеграл:

In [30]:
quad(lambda x,y:1/(1+x*y),[0,1],[0,1])
Out[30]:
mpf('0.82246703342411321823620758332301259460947495060339899')

Сумма ряда

In [31]:
with extradps(5):
    s=nsum(lambda n:(-1)**(n-1)/n**4,(1,inf))
    print(s)

0.9470328294972459175765032344735219149279070829288860442
In [32]:
pslq([s,pi**4])
Out[32]:
[-720, 7]

То есть эта сумма, вероятно, равна $\frac{7}{720} \pi^4$.

Дифференциальные уравнения

In [33]:
a=mpf('0.2')
def f(t,x):
    global a
    return [x[1],-x[0]-2*a*x[1]]
In [34]:
x=odefun(f,0,[1,0])
In [35]:
x(1)
Out[35]:
[mpf('0.59496623263788777500734762840237378880987158574261147'),
 mpf('-0.69387986210972080683187214187798497495035321299936026')]
In [36]:
plot([lambda t:x(t)[0],lambda t:x(t)[1]],[0,10])

Матрицы

Матрицы разреженные, реализованы как словари. Квадратная матрица

In [37]:
matrix(2)
Out[37]:
matrix(
[['0.0', '0.0'],
 ['0.0', '0.0']])

Прямоугольная матрица

In [38]:
M=matrix(2,3)
M
Out[38]:
matrix(
[['0.0', '0.0', '0.0'],
 ['0.0', '0.0', '0.0']])
In [39]:
M.rows,M.cols
Out[39]:
(2, 3)
In [40]:
M[0,1]=1
M
Out[40]:
matrix(
[['0.0', '1.0', '0.0'],
 ['0.0', '0.0', '0.0']])

Операции с матрицами

In [41]:
M1=matrix([[0,1],[1,0]])
M2=matrix([[1,0],[0,-1]])
In [42]:
M1+M2
Out[42]:
matrix(
[['1.0', '1.0'],
 ['1.0', '-1.0']])
In [43]:
M1*M2
Out[43]:
matrix(
[['0.0', '-1.0'],
 ['1.0', '0.0']])
In [44]:
M2*M1
Out[44]:
matrix(
[['0.0', '1.0'],
 ['-1.0', '0.0']])
In [45]:
M1**(-1)
Out[45]:
matrix(
[['0.0', '1.0'],
 ['1.0', '0.0']])

Решение системы линейных уравнений

In [46]:
A=matrix([[1,2],[3,4]])
b=matrix([1,-1])
b
Out[46]:
matrix(
[['1.0'],
 ['-1.0']])
In [47]:
x=lu_solve(A,b)
x
Out[47]:
matrix(
[['-3.0'],
 ['2.0']])
In [48]:
A*x-b
Out[48]:
matrix(
[['0.0'],
 ['0.0']])

Собственные значения и собственные векторы

In [49]:
l,u=eig(A)
l
Out[49]:
[mpf('-0.37228132326901432992530573410946465911013222899139797'),
 mpf('5.3722813232690143299253057341094646591101322289914067')]
In [50]:
u
Out[50]:
matrix(
[['-0.82456484013239376536905071707877267896095335074304', '-0.42222915041526045335929057758178658089159736117701'],
 ['0.56576746496899228472288762798052673125191630934726', '-0.92305231425019333318861560854941073593095247730112']])

Диагональная матрица

In [51]:
L=diag(l)
L
Out[51]:
matrix(
[['-0.3722813232690143299253057341094646591101322289914', '0.0'],
 ['0.0', '5.3722813232690143299253057341094646591101322289914']])
In [52]:
A*u-u*L
Out[52]:
matrix(
[['-2.0045735325691467346054023503636114091113976600788e-51', '5.3455294201843912922810729343029637576303937602101e-51'],
 ['0.0', '1.069105884036878258456214586860592751526078752042e-50']])