时间戳

首页 > 解决方案 > 自动将算术运算四舍五入到小数点后八位

python-3.x - 自动将算术运算四舍五入到小数点后八位

问题描述

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

或者，是否有可以指定算术精度的语言/程序？

标签： python-3.xnumpynumerical-methodsnumerical-analysis

解决方案

计算期间的精度/舍入不太可能成为问题。

为了测试这一点，我使用包含您目标精度的精度运行计算：一次使用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在低精度下没有表现出我的预期，所以我放弃了它。

推荐阅读