Блог

Я выполняю некоторое упражнение по числовому анализу, в котором мне нужно рассчитать решение линейной системы с использованием определенного алгоритма. Мой ответ отличается от ответа книги на несколько десятичных знаков, что, как я полагаю, связано с ошибками округления. Есть ли способ, с помощью которого я могу автоматически устанавливать арифметику на округление до восьми знаков после запятой после каждой арифметической операции? Ниже приведен мой код на python.

 import numpy as np
A1 = [4, -1, 0, 0, -1, 4, -1, 0,
     0, -1, 4, -1, 0, 0, -1, 4]
A1 = np.array(A1).reshape([4,4])
I = -np.identity(4)
O = np.zeros([4,4])

A = np.block([[A1, I, O, O],
             [I, A1, I, O],
             [O, I, A1, I],
             [O, O, I, A1]])

b = np.array([1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6])

def conj_solve(A, b, pre=False):
    n = len(A)
    
    C = np.identity(n)        
    
    if pre == True:
        for i in range(n):
            C[i, i] = np.sqrt(A[i, i])
            
    Ci = np.linalg.inv(C)
    Ct = np.transpose(Ci)
    
    x = np.zeros(n)
    r = b - np.matmul(A, x)
    w = np.matmul(Ci, r)
    v = np.matmul(Ct, w)
    alpha = np.dot(w, w)
    
    for i in range(MAX_ITER):
        if np.linalg.norm(v, np.infty) < TOL:            
            print(i 1, "steps")            
            print(x)
            print(r)
            return
        u = np.matmul(A, v)
        t = alpha/np.dot(v, u)
        x = x   t*v
        r = r - t*u
        w = np.matmul(Ci, r)
        beta = np.dot(w, w)
        
        if np.abs(beta) < TOL:
            if np.linalg.norm(r, np.infty) < TOL:
                print(i 1, "steps")
                print(x)
                print(r)
                return
        s = beta/alpha
        v = np.matmul(Ct, w)   s*v
        alpha = beta
    print("Max iteration exceeded")
    return x

MAX_ITER = 1000
TOL = 0.05

sol = conj_solve(A, b, pre=True)
  

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

ИЛИ, есть ли язык / программа, где я могу указать точность арифметики?

Точность / округление во время вычисления вряд ли будет проблемой.

Чтобы проверить это, я выполнил вычисление с точностью, заключающей в скобки точность, к которой вы стремитесь: один раз с np.float64 , и один раз с np.float32 . Вот таблица напечатанных результатов, их приблизительная десятичная точность и результат вычисления (т.Е. Первое напечатанное значение массива).

 numpy type       decimal places         result
-------------------------------------------------
np.float64               15              2.55516527
np.float32                6              2.5551653
  

Учитывая, что они настолько согласуются, я сомневаюсь, что промежуточная точность в 8 знаков после запятой даст ответ, который не находится между этими двумя результатами (т. 2.55613420 Е. отключен в 4-й цифре).

Это не часть, не часть моего ответа, а комментарий к использованию mpmath . Спрашивающий предложил это в комментариях, и это была и моя первая мысль, поэтому я провел быстрый тест, чтобы проверить, работает ли он так, как я ожидал, при вычислениях с низкой точностью. Этого не произошло, поэтому я отказался от него (но я не эксперт в этом).

Вот моя тестовая функция, в основном умножающаяся 1/N на N и 1/N повторно, чтобы подчеркнуть ошибку 1/N .

 def precision_test(dps=100, N=19, t=mpmath.mpf):
    with mpmath.workdps(dps):  
        x = t(1)/t(N)
        print(x)
        y = x
        for i in range(10000):
            y *= x
            y *= N
        print(y)
  

Это работает, как и ожидалось, например, np.float32 :

 precision_test(dps=2, N=3, t=np.float32)
# 0.33333334
# 0.3334327041164994
  

Обратите внимание, что ошибка распространилась на более значащие цифры, как и ожидалось.

Но mpmath я никогда не мог добиться этого (тестирование с диапазоном dps и различными простыми N значениями):

 precision_test(dps=2, N=3)
# 0.33
# 0.33
  

Из-за этого теста я решил mpmath , что он не даст нормальных результатов для вычислений с низкой точностью.

TL; DR:
mpmath не вел себя так, как я ожидал, с низкой точностью, поэтому я отказался от него.