问题描述
我正在做一些数值分析练习，我需要使用specific algorithm计算线性系统的解。我的答案与本书的答案有些不同，我认为这是四舍五入引起的。有没有一种方法可以在每次算术运算后自动将算术设置为四舍五入到小数点后八位？以下是我的python代码。
import numpy as np
A1 = [4,-1,4,\
0,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,[O,I],A1]])
b = np.array([1,2,3,5,6,7,8,9,1,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,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位）。 / p>
这不是我的答案的一部分，但这是对使用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的表现不符合我期望的那样，因此我放弃了它。