曲率数值计算：数值结果与解析值偏差问题及解决方法

首先要指出一个关键误区：你提到的解析解公式-0.02*(x-500)**2 + 250其实是原曲线的方程，而不是曲率的解析解——这正是二者看起来偏差极大的核心原因！我们先推导正确的曲率解析解，再修正你的数值计算代码。

第一步：推导正确的曲率解析解 #

你的原曲线方程是：

y = -0.02*(x-500)**2 + 250

对x求导得到：

• 一阶导数：dy/dx = -0.04*(x-500)
• 二阶导数：d²y/dx² = -0.04

代入平面曲线曲率公式：

curvature = |d²y/dx²| / (1 + (dy/dx)²)^(3/2)

代入导数后得到解析曲率：

cur_analytical = 0.04 / (1 + (0.04*(x-500))**2)**1.5

这才是你应该用来对比的基准值。

第二步：修正数值计算代码的问题 #

你的代码有两个主要问题，导致数值曲率和解析解偏差：

1. 导数计算未指定x轴步长 #

np.gradient默认步长为1，但如果你的x轴实际间隔不是1，直接计算dy = np.gradient(data[:,1])会得到错误的一阶导数数值，进而影响后续所有计算。必须显式传入x轴的均匀步长。

2. 二阶导数的边缘精度（可选优化） #

直接嵌套np.gradient计算二阶导数时，边缘点会使用一阶差分近似，精度略低。对于均匀间隔数据，可以用二阶中心差分提升中间点的精度。

修正后的完整代码 #

import numpy as np

# 加载数据
data = np.loadtxt('newsorted.txt')
x = data[:, 0]
y = data[:, 1]

# 获取x轴的均匀步长（因为是等间距数据，直接取相邻两点差）
h = x[1] - x[0]

# 计算一阶导数 dy/dx：显式指定步长h
dy = np.gradient(y, h)

# 计算二阶导数 d²y/dx²
# 方法1：用np.gradient嵌套（简单易用，边缘精度一般）
d2y = np.gradient(dy, h)

# 方法2：二阶中心差分（中间点精度更高，边缘点可单独处理）
# d2y = np.zeros_like(y)
# d2y[1:-1] = (y[2:] - 2*y[1:-1] + y[:-2]) / (h**2)
# # 边缘点用前向/后向差分的二阶近似
# d2y[0] = (y[2] - 2*y[1] + y[0]) / (h**2)
# d2y[-1] = (y[-1] - 2*y[-2] + y[-3]) / (h**2)

# 计算曲率
cur = np.abs(d2y) / (1 + dy**2)**1.5

第三步：验证对比 #

用下面的代码计算解析曲率，再和数值结果对比：

# 计算解析曲率
dy_analytical = -0.04*(x - 500)
cur_analytical = 0.04 / (1 + dy_analytical**2)**1.5

此时你会发现，数值计算的曲率和解析解几乎完全重合，边缘点的微小误差可以通过优化边缘差分方法进一步减小。

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