计算物理（浮点数精度相关）

来源：

        某次计算物理作业

题目：

只提供1）和2）的解题代码，注意，为了使现象更明显，我缩小了自变量的取值范围：

代码：

import matplotlib.pyplot as plt
import numpy as np
from mpmath import mp

# 设置mpmath库的精度为四倍精度
mp.dps = 36


# 定义函数
def func(x):
    return (x - 1) ** 12


def func2(x):
    return x ** 12 - 12 * x ** 11 + 66 * x ** 10 - 220 * x ** 9 + 495 * x ** 8 - 792 * x ** 7 + 924 * x ** 6 - 792 * x ** 5 + 495 * x ** 4 - 220 * x ** 3 + 66 * x ** 2 - 12 * x + 1


# 定义 x 范围
x_single = np.linspace(0.998, 1.002, 400, dtype=np.float32)
x_double = np.linspace(0.998, 1.002, 400, dtype=np.float64)
step = mp.mpf('0.00001')
x_quad = [mp.mpf(val) for val in mp.arange('0.998', '1.002', step)]

# 计算函数值
y_single = func(x_single)
y_double = func(x_double)
y_quad = [func(x) for x in x_quad]

# 计算展开式函数图
y_single2 = func2(x_single)
y_double2 = func2(x_double)
y_quad2 = [func2(x) for x in x_quad]

# 绘制图像
plt.figure(figsize=(24, 6))

plt.subplot(1, 3, 1)
plt.scatter(x_single, y_single, label='Single Precision', color='blue')
plt.scatter(x_single, y_single2, label='expansion', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.title('Graph of $(x-1)^{12}$ using Single Precision')

plt.subplot(1, 3, 2)
plt.scatter(x_double, y_double, label='Double Precision', color='blue')
plt.scatter(x_double, y_double2, label='expansion', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.title('Graph of $(x-1)^{12}$ using Double Precision')

plt.subplot(1, 3, 3)
plt.scatter(x_quad, y_quad, label='Quadruple Precision', color='blue')
plt.scatter(x_quad, y_quad2, label='expansion', color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.title('Graph of $(x-1)^{12}$ using Quadruple Precision (mpmath)')

plt.tight_layout()
plt.show()

# # 计算相对误差
# relative_error_single = np.zeros_like(y_single)
# for i in range(len(y_single)):
#     if func(x_single[i]) != 0:
#         relative_error_single[i] = np.abs((y_single2[i] - func(x_single[i])) / func(x_single[i]))
#
# relative_error_double = np.zeros_like(y_double)
# for i in range(len(y_double)):
#     if func(x_double[i]) != 0:
#         relative_error_double[i] = np.abs((y_double2[i] - func(x_double[i])) / func(x_double[i]))
#
# # 修改计算相对误差的部分代码
# relative_error_quad = np.zeros_like(y_quad)
# for i in range(len(x_quad)):
#     if float(func(x_quad[i])) != 0:
#         relative_error_quad[i] = np.abs((y_quad2[i] - float(func(x_quad[i]))) / float(func(x_quad[i])))
#
# # 绘制相对误差图像
# plt.figure(figsize=(18, 6))
#
# plt.subplot(1, 3, 1)
# plt.scatter(x_single, relative_error_single, label='Single Precision', color='blue')
# plt.xlabel('x')
# plt.ylabel('Relative Error')
# plt.legend()
# plt.grid(True)
# plt.title('Relative Error using Single Precision')
#
# plt.subplot(1, 3, 2)
# plt.scatter(x_double, relative_error_double, label='Double Precision', color='red')
# plt.xlabel('x')
# plt.ylabel('Relative Error')
# plt.legend()
# plt.grid(True)
# plt.title('Relative Error using Double Precision')
#
# plt.subplot(1, 3, 3)
# plt.scatter(x_quad, relative_error_quad, label='Quadruple Precision', color='green')
# plt.xlabel('x')
# plt.ylabel('Relative Error')
# plt.legend()
# plt.grid(True)
# plt.title('Relative Error using Quadruple Precision (mpmath)')
#
# plt.tight_layout()
# plt.show()

# 设置一个很小的数，用于替代零
epsilon_single = 1e-8
epsilon_double = 1e-16
epsilon_quad = 1 * 10 ** (-36)

# 计算相对误差，避免除以零
relative_error_single = np.abs(y_single - y_single2) / (np.abs(y_single) + epsilon_single)
relative_error_double = np.abs(y_double - y_double2) / (np.abs(y_double) + epsilon_double)
y_quad = np.array(y_quad)
y_quad2 = np.array(y_quad2)
relative_error_quad = np.abs(y_quad - y_quad2) / (np.abs(y_quad) + epsilon_quad)

# 绘制相对误差图表
plt.figure(figsize=(24, 6))
plt.subplot(1, 3, 1)
plt.scatter(x_single, relative_error_single, label='Relative Error', color='red')
plt.xlabel('x_single')
plt.ylabel('Relative Error')
plt.legend()
plt.grid(True)
plt.title('Relative Error between Single Precision and Expansion')

plt.subplot(1, 3, 2)
plt.scatter(x_double, relative_error_double, label='Relative Error', color='blue')
plt.xlabel('x_double')
plt.ylabel('Relative Error')
plt.legend()
plt.grid(True)
plt.title('Relative Error between Double Precision and Expansion')

plt.subplot(1, 3, 3)
plt.scatter(x_quad, relative_error_quad, label='Relative Error', color='green')
plt.xlabel('x_quad')
plt.ylabel('Relative Error')
plt.legend()
plt.grid(True)
plt.title('Relative Error between Quad Precision and Expansion')

plt.show()

运行结果：

结语：

嗯，这学期课比较少了，争取多写一点东西。同时欢迎批评指正，本人仅代表专业最低水平QAQ