在Matlab中用矩阵微分算子计算离散一阶导数及异常排查

振幅缩放错误的原因 #

你犯了一个基础的导数公式错误：你的目标函数是f = sin(k.*x)，根据链式法则，它的导数应该是k*cos(k.*x)，而非你以为的cos(x)。代入你的参数计算，k = 2πf/c ≈ 18.326，所以导数的振幅是k，不是1。你当前用中心差分得到的结果其实是在近似k*cos(kx)，只是因为你错误设定了预期值，才会觉得振幅存在缩放问题。

你可以在代码中添加理论导数的绘制来验证：

plot(k*cos(k.*x), '--r')

对比后会发现差分结果的振幅和理论值完全匹配。

首尾样本偏离预期的原因 #

你构建的微分矩阵D采用的是中心差分格式，这个格式的核心公式是f’(x_i) ≈ (f(x_{i+1}) - f(x_{i-1}))/(2dx)，但它只适用于序列的中间点（第2到第N-1个样本）。

对于首尾两个样本：

• 第一个样本没有前一个点x_0，无法使用中心差分，但你的D矩阵第一行是[0, 1/(2dx), 0, ..., 0]，计算的是f(2)/(2dx)，这完全不符合导数的近似逻辑；
• 最后一个样本没有后一个点x_{N+1}，你的D矩阵最后一行是[0, ..., 0, -1/(2dx), 0]，计算的是-f(N-1)/(2dx)，同样是错误的近似。

修正方案 #

对首尾点分别使用前向差分和后向差分，中间点保留中心差分，修改后的矩阵构造代码如下：

N = length(f);
D = zeros(N, N);
% 中间点：中心差分
for i = 2:N-1
    D(i, i-1) = -1/(2*dx);
    D(i, i+1) = 1/(2*dx);
end
% 首点：前向差分
D(1, 1) = -1/dx;
D(1, 2) = 1/dx;
% 末点：后向差分
D(N, N-1) = -1/dx;
D(N, N) = 1/dx;

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