Умножение матриц дает необычный результат в Python (SciPy / PyLab) - PullRequest
Новая
       35

Умножение матриц дает необычный результат в Python (SciPy / PyLab)

5 голосов
/ 06 февраля 2011

Я новичок в Python, и немного устала от своей линейной алгебры, так что, возможно, это простой вопрос. Я пытаюсь реализовать расширение серии Тейлора на матрице для вычисления exp (A), где A - простая матрица 3x3. Формула, кстати, для этого расширения является сумма (A ^ n / n!).

Моя процедура работает нормально до n = 9, но при n = 10 числа в Матрице внезапно становятся отрицательными. Это проблема.

A ** 9 матрица ([[250130371, 506767656, 688136342], [159014912, 322268681, 437167840], [382552652, 775012944, 1052574077]])

A ** 10 матрица ([[- 1655028929, 1053671123, -1327424345], [1677887954, -895075635, 319718665], [-257240602, -409489685, -1776533068]])

Интуитивно, A ^ 9 * A должен производить большие числа для каждого члена матрицы, но, как вы можете видеть, A ^ 10 не дает такого результата.

Есть идеи? 

from scipy import *
from numpy import *
from numpy.linalg import *
#the matrix I will use to implement exp(A)
A = mat('[1 3 5; 2 5 1; 2 3 8]')
#identity matrix
I = mat('[1 0 0; 0 1 0; 0 0 1]')
#first step in Taylor Expansion (n=0)
B = I
#second step in Taylor Expansion (n=1)
B += A
#start the while loop in the 2nd step
n = 2
x=0
while x<10:
    C = (A**n)/factorial(n)
    print C
    print " "
    n+=1
    B+= C
    print B
    x+=1

print B

Спасибо за любую помощь, которую вы можете оказать!

Ответы [ 2 ]

8 голосов
/ 06 февраля 2011

Ваша матрица создана с элементами типа int32 (32-разрядное целое число). Вы можете увидеть это, напечатав значение A.dtype. 32-разрядные целые числа могут содержать значения только до 2 миллиардов, поэтому после этого они обернутся в отрицательные значения.

Если 64-разрядные целые числа достаточно велики, вы можете использовать их вместо: 

A = mat('[1 3 5; 2 5 1; 2 3 8]', dtype=numpy.int64)

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

A = mat('[1 3 5; 2 5 1; 2 3 8]', dtype=float)

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

3 голосов
/ 06 февраля 2011

Я не знаю много о научном питоне, но я знаю, что происходит не так. Кажется, что элементы матрицы представлены как 32-битные целые числа со знаком. Это означает, что они ограничены диапазоном -2 ^ 31 <= x <2 ^ 31. На A ^ 10 числа становятся слишком большими, и они «оборачиваются». Я проверил, и верхний левый коэффициент на самом деле равен 2639938267, что при обтекании дает 2639938267 - 2 ^ 32 = -1655028929. </p>

Я не знаю, как установить тип данных в Python, поэтому я не знаю, как вы можете решить эту проблему. Я уверен, что это возможно, хотя.

(Я также мог бы предложить вам попробовать sage: www.sagemath.org, математическое программное обеспечение, основанное на python. Оно автоматически использует бесконечную точность. Именно так я и проверил эти числа только сейчас).

Удачи! Timo

...