Matplotlib Quiver绘制引力矢量场时近质点位置异常问题

我尝试用Matplotlib的quiver函数绘制引力矢量场，编写了一个可计算任意数量质点所产生引力场的函数。该函数可为矢量绘图返回场的x、y分量以及对应的x、y坐标，还有一个后续代码会用到的二维场数组（绘图无需该数组）。x、y分量等四个列表的对应索引元素为同一场点的数据，即索引0处可获取该点的场分量及位置。

当前代码如下：

import matplotlib.pyplot as plt
import numpy as np
import math

G = 6.674 * 10**(-11)

# 计算所有质点在屏幕各点产生的引力场
def gravfield(mass, x, y, screen_width, screen_height, give_components = False):

    # 定义屏幕坐标点
    x_screen = list(range(screen_width))
    y_screen = list(range(screen_height))

    # 初始化场强数组、场分量列表及坐标列表
    field = np.zeros((screen_width, screen_height))
    field_x = []
    field_y = []
    field_x_position = []
    field_y_position = []
        
    # 用矢量法计算每个质点产生的引力场叠加
    for i in range(len(x_screen)):
        for n in range(len(y_screen)):

            # 初始化当前点的合场强分量
            x_field = 0
            y_field = 0
            
            for j in range(len(mass)):
                
                # 计算质点到场点的位置矢量
                x_distance = x[j] - x_screen[i]
                y_distance = y[j] - y_screen[n]
                distance = np.sqrt((x_distance**2) + (y_distance**2))
                
                if not math.isclose(distance, 0):
                        # 计算单位矢量
                        x_unit_vector = x_distance / distance
                        y_unit_vector = y_distance / distance

                        # 计算场强大小
                        magnitude = (G * mass[j]) / (distance ** 2)

                        # 叠加各质点的场强分量
                        x_field += magnitude * x_unit_vector
                        y_field += magnitude * y_unit_vector
                else:
                    pass
                
        # 将当前点的场强数据存入对应结构
            if give_components == True:
                field_x.append(x_field)
                field_y.append(y_field)
                field[n, i] = np.sqrt(x_field**2 + y_field**2)
                field_x_position.append(i)
                field_y_position.append(n)
            else:
                field[n, i] = np.sqrt(x_field**2 + y_field**2)

    # 返回指定结果
    if give_components == True:
        return field, field_x, field_y, field_x_position, field_y_position
    else:
        return field

# 测试参数
mass = [4e23, 8e23]
x_cord = [10, 15]
y_cord = [10, 15]
screen_width = 20
screen_height = 20
give_components = True
field, field_x, field_y, field_x_position, field_y_position = gravfield(mass, x_cord, y_cord, screen_width, screen_height, give_components)

# 绘制矢量场
fig, ax = plt.subplots()
ax.quiver(field_x_position, field_y_position, field_x, field_y)

# 设置坐标轴范围
ax.set_xlim([min(field_x_position), max(field_x_position)])
ax.set_ylim([min(field_y_position), max(field_y_position)])

plt.show()

该函数用于场强的颜色编码显示时正常，但使用quiver函数绘制时出现异常：除质点附近区域外，场表现正常，而近质点处出现了方向混乱、长度明显异于周围的矢量。已确认不是除零错误，计算时已忽略距离为零的情况。

异常绘图表现：两个质点（坐标(10,10)和(15,15)）附近的矢量方向杂乱，长度远大于或小于周边正常矢量，其余区域的矢量均指向对应质点，符合引力场的预期分布。

内容的提问来源于stack exchange，提问作者Valter Somlai