我试图重复我在一篇论文中找到的一个例子。

我必须解决这个 ODE：

$25 \ddot{x} + 15 \dot{x} + 330000 x = p(t)$

其中是频带限制在 10-25 Hz 范围内的白噪声序列。$p(t)$

然后它说系统是使用 Runge-Kutta 程序模拟的。采样频率设置为 1000 Hz，并将高斯白噪声添加到数据中，使噪声占信号 rms 值的 5%（如何使用最后的信息？）。

主要问题与带限白噪声有关。我按照我在这里找到的说明https://dsp.stackexchange.com/questions/9654/how-to-generate-band-limited-gaussian-white-noise并编写了以下代码：

% White noise excitation
% design FIR filter to filter noise in the range [10-25] Hz
Fs = 1000; % sampling frequency

% Time infos (for white noise)
T = 1/Fs;
tspan = 0:T:4; % I saw from the plots in the paper that the simulation last 4 seconds
tstart = tspan(1);
tend = tspan (end);

b = fir1(64, [10/(Fs/2) 25/(Fs/2)]);
% generate Gaussian (normally-distributed) white noise
noise = randn(length(tspan), 1);
% apply filter to yield bandlimited noise
p = filter(b,1,noise);

现在我必须为 ode 定义函数，但我不知道如何给它（白噪声）：$p$

我试过这样编码

odefun = @(t,u,m,k1,c,p)[u(2); 1/m*(-k1*u(1)-c*u(2)+p)];
q0 = [0;0]; % Initial conditions
[t,q,te,qe,ie] = ode45(@(t,u)odefun(t,u,m,k1,c,p),tspan,q0(:,1),options);

但我收到此错误消息：

@(T,U)ODEFUN(T,U,M,K1,C,NB) 返回一个长度为 4002 的向量，但初始条件向量的长度为 2。@(T,U)ODEFUN(T 返回的向量,U,M,K1,C,NB) 和初始条件向量必须具有相同数量的元素。

我已经看过以下问题，但我不明白问题可能出在哪里。

https://stackoverflow.com/questions/14343563/how-solve-a-system-of-ordinary-differntial-equation-with-time-dependent-paramete?lq=1

https://stackoverflow.com/questions/19732442/solving-an-ode-when-the-function-is-given-as-discrete-values-matlab?lq=1

https://stackoverflow.com/questions/19732442/solving-an-ode-when-the-function-is-given-as-discrete-values-matlab

编辑

我可能已经找到了一种方法：

function [y] = p(t)
    b = fir1(64, [10/(Fs/2) 25/(Fs/2)]);
    % generate Gaussian (normally-distributed) white noise
    noise = randn(length(t), 1);
    % apply filter to yield bandlimited noise
    y = filter(b,1,noise);
end

odefun = @(t,u,m,k1,c,p)[u(2); 1/m*(-k1*u(1)-c*u(2)+p(t))];

[t,q,te,qe,ie] = ode45(@(t,u)odefun2(t,u,m,k2,c,@p),tspan,q0(:,1),options);

你觉得好看吗？