原始结果
运行您提供的代码并将结果与your other question 中发布的 RK4 代码进行比较，我们得到以下响应：
蓝色曲线代表基于 FFT 的实现，红色曲线代表您的替代 RK4 实现。正如您所指出的，曲线完全不同。
得到正确的响应
最明显的问题当然是幅度，而本题贴出的代码幅度差异的主要来源与我在my answer to your other question中指出的相同：
- RK4 实现执行数字积分，通过积分步骤
dt 正确缩放总和值。基于 FFT 的实现缺乏这种缩放。
- 在基于 FFT 的实现中使用的脉冲响应与按质量
m 缩放的驱动力一致，这是 RK4 实现中缺少的一个因素。
解决这两个问题会产生更接近但仍不相同的响应。正如您可能发现的那样，鉴于您发布的其他问题的代码发生了变化，缺少的另一件事是输入和脉冲响应的零填充，没有它，您将获得 圆形 卷积而不是比线性卷积：
f=[cos(t),zeros(1,length(t)-1)]; %force f
h=[1/m/wd*exp(-z*wn*t).*sin(wd*t),zeros(1,length(t)-1)]; %impulse response
最后，确保卷积产生良好结果的最后一个元素是使用对无限长脉冲响应的良好近似。多长时间足够长取决于脉冲响应的衰减率。使用您提供的参数，在大约11*pi 之后，脉冲响应将下降到其原始值的 1%。因此，将时间跨度扩展到tMax=14*pi（在脉冲响应消失后包括完整的2*pi 循环）可能就足够了。
获得稳态响应
获得稳态响应的最简单方法是丢弃初始瞬态。在这个过程中，我们丢弃了参考驱动力的整数个周期（这当然需要知道驱动力的基频）：
T0 = tMax-2*pi;
delay = min(find(t>T0));
sol = sol(delay:end);
plot([0:length(sol)-1]*dt, sol, 'b');
axis([0 2*pi]);
然后得到的响应是：
蓝色曲线再次代表基于 FFT 的实现，红色曲线代表您的替代 RK4 实现。好多了！
另一种方法
计算多个周期的响应，等待瞬态响应消退，然后提取对应的剩余样本
尽管由于 FFT 计算仍然相当快，但到稳定状态可能看起来有点浪费。
所以，让我们回过头来看看问题域。你可能知道，
mass-spring-damper system 由微分方程控制：
f(t) 在这种情况下是驱动力。
请注意，齐次方程的通解具有以下形式：
然后关键是要认识到c>0 和m>0 在稳定状态下消失的情况下的一般解决方案（t 趋于无穷大）。
因此，稳态解仅取决于非齐次方程的特定解。
这个特殊的解决方案可以通过method of undetermined coefficients找到，求驱动的形式
相应地假设解决方案具有以下形式
代入微分方程得到方程：
因此，解决方案可以实现为：
EF0 = [wn*wn-w*w 2*z*wn*w; -2*z*wn*w wn*wn-w*w]\[1/m; 0];
sol = EF0(1)*cos(w*t)+EF0(2)*sin(w*t);
plot(t, sol);
在你的情况下w=2*pi。
概括
通过将驱动力表示为
傅里叶级数（假设驱动力函数满足Dirichlet conditions）：
可以相应地假设特解具有以下形式
求解特定解的方法与前面的情况非常相似。这导致以下实现：
% normalize
y = F/m;
% compute coefficients proportional to the Fourier series coefficients
Yw = fft(y);
% setup the equations to solve the particular solution of the differential equation
% by the method of undetermined coefficients
k = [0:N/2];
w = 2*pi*k/T;
A = wn*wn-w.*w;
B = 2*z*wn*w;
% solve the equation [A B;-B A][real(Xw); imag(Xw)] = [real(Yw); imag(Yw)] equation
% Note that solution can be obtained by writing [A B;-B A] as a scaling + rotation
% of a 2D vector, which we solve using complex number algebra
C = sqrt(A.*A+B.*B);
theta = acos(A./C);
Ywp = exp(j*theta)./C.*Yw([1:N/2+1]);
% build a hermitian-symmetric spectrum
Xw = [Ywp conj(fliplr(Ywp(2:end-1)))];
% bring back to time-domain (function synthesis from Fourier Series coefficients)
x = ifft(Xw);
最后一点
我在上面的推导中故意避免了无阻尼的c=0 情况。在这种情况下，振荡永远不会消失，齐次方程的通解也不一定是微不足道的。
在这种情况下，最终的“稳态”可能与驱动力具有相同的周期，也可能不同。事实上，如果通解的周期振荡与驱动力的周期通过有理数（整数比）无关，则它可能根本不是周期性的。