自动将算术运算四舍五入到八位小数

Question

我正在做一些数值分析练习，我需要使用 specific algorithm 计算线性系统的解。我的答案与书中的答案有一些小数点的不同，我认为这是由于舍入误差造成的。有没有一种方法可以让我在每次算术运算后自动将算术设置为四舍五入小数点后八位？以下是我的 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。

或者，是否有 language/program 可以指定算术精度的地方？

Answer 1

Precision/rounding 在计算过程中不太可能是问题。

为了测试这一点，我运行计算的精度涵盖了您的目标精度：一次使用 np.float64，一次使用 np.float32。这是打印结果的table，它们的近似小数精度，以及计算结果（即第一个打印的数组值）。

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 的运行ge 和各种素数 N 值进行测试）：

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

因为这个测试，我决定 mpmath 不会为低精度计算提供正常结果。

TL;DR:
mpmath 在低精度下没有达到我的预期，所以我放弃了它。

自动将算术运算四舍五入到八位小数

Automatically round arithmetic operations to eight decimals

python-3.x

numpy

numerical-methods

numerical-analysis