时间戳

首页 > 解决方案 > 在笛卡尔网格上计算的函数的数值径向导数

python - 在笛卡尔网格上计算的函数的数值径向导数

问题描述

我有一个在 3D 笛卡尔网格上评估的径向对称函数。如何数值计算函数的径向导数？

对于一个简单的例子（球面高斯），计算导数 df/dx、df/dy 和 df/dz：

# Parameters
start = 0
end = 5
n = 20

# Variables
x = np.linspace(start, end, num=n)
y = np.linspace(start, end, num=n)
z = np.linspace(start, end, num=n)
dx = (end - start) / n
dy = (end - start) / n
dz = (end - start) / n
x_grid, y_grid, z_grid = np.meshgrid(x, y, z)
eval_xyz = np.exp(-(x_grid ** 2 + y_grid ** 2 + z_grid ** 2))

# Allocate
df_dx = np.zeros((n, n, n))
df_dy = np.zeros((n, n, n))
df_dz = np.zeros((n, n, n))

# Calculate Cartesian gradient numerically
for x in range(eval_xyz.shape[0] - 1):
    for y in range(eval_xyz.shape[1] - 1):
        for z in range(eval_xyz.shape[2] - 1):

            df_dx[x, y, z] = (eval_xyz[x + 1, y, z] - eval_xyz[x, y, z]) / dx
            df_dy[x, y, z] = (eval_xyz[x, y + 1, z] - eval_xyz[x, y, z]) / dy
            df_dz[x, y, z] = (eval_xyz[x, y, z + 1] - eval_xyz[x, y, z]) / dz

那么是否可以从笛卡尔导数轻松计算径向导数 df/dr？

标签： pythonnumpygeometry

解决方案

诀窍是将径向导数表示为笛卡尔导数的总和，同时考虑每个点的 theta 和 phi，这可以用笛卡尔坐标表示为：方程

因此代码变为：

theta_val = theta(i * dx, j * dy, k * dz)
phi_val = phi(i * dx, j * dy)
            
df_dr[i, j, k] = df_dx[i, j, k] * np.sin(theta_val) * np.cos(phi_val) \
                 + df_dy[i, j, k] * np.sin(theta_val) * np.sin(phi_val) \
                 + df_dz[i, j, k] * np.cos(theta_val)

仔细计算 theta 和 phi 以处理除以零

def theta(x, y, z):

    if x == 0 and y == 0 and z == 0:
        return 0

    elif z == 0:
        return np.pi / 2

    elif x == 0 and y == 0:
        return 0

    else:
        return np.arctan(np.sqrt(x ** 2 + y ** 2) / z)


def phi(x, y):

    if x == 0 and y == 0:
        return 0

    elif x == 0:
        return np.pi / 2

    elif y == 0:
        return 0

    else:
        return math.atan2(y, x)

推荐阅读