MATLAB二阶微分方程离散解绘制功率谱

在MATLAB中写出图中的二阶微分方程，通过四阶龙格库塔法求解x(t),y(t)的离散解，离散解结果无误，但是根据离散解通过写的傅里叶变换程序求功率谱如左图，始终无法画出右图中的结果。而且画出来的结果一直受计算时长的影响较大，但参考文献中没有写出。想请问应该写怎么样的程序才能得出右图的功率谱。

img

clear;
clc
fwind=0.01144;
R1=3.17;
R3=-1.61;
R5=-2.76;
g=9.81;
I=5.37*10^10;    
delta=1.47*10^8; 
r1=delta*g*R1/I;
r3=delta*g*R3/I;
r5=delta*g*R5/I;

%options=odeset('maxstep',0.01);%%时间步长调整
%[t,y]=ode45(@vdp1,[-30 30],[0.875;0],options);%%四阶龙哥库塔法求解微分方程组，y(:,1)即x的解，y(:,2)即y的解
odefun2 = @(t,y)[y(2);-(r1*y(1)+r3*y(1)^3+r5*y(1)^5)+fwind*(1-y(1)^2/2+y(1)^4/24)]; % 化二阶常微分方程为二阶常微分方程组

ts = -30; %计算开始时间
te = 30; %计算结束时间
h = 0.001;  %时间步长
t = ts:h:te; %计算数据点
y = zeros(length(t),2); 
y(1,:)=[0.875,0]; %参数初值

for n = 1:length(t)-1
k1 = odefun2(t(n),y(n,:)');
k2 = odefun2(t(n)+h/2,y(n,:)'+h/2*k1); %这里要特别注意同型向量相加
k3 = odefun2(t(n)+h/2,y(n,:)'+h/2*k2);
k4 = odefun2(t(n)+h,y(n,:)'+h*k3);
y(n+1,:) = y(n,:)+h/6*(k1+2*k2+2*k3+k4)';%四级四阶Runge-Kutta公式
end
%%傅里叶变换
N=length(y(:,2));
y_OMEGA=fft(y(:,2));
f=(0:N-1)/(N*h);%%采样频率

omega_0=(delta*g*R1/I)^0.5;  
omega=f*2*pi;
OMEGA=omega/omega_0; %求出相应横坐标Ω

S_y=(abs(y_OMEGA)).^2./(N); %Sy(Ω)
plot(OMEGA,S_y);
axis([0 1.5, -inf,inf]);

麻烦复制一下代码

参考GPT和自己的思路，在 MATLAB 中，可以通过以下方式编写二阶微分方程，使用四阶龙格-库塔法求解离散解：

% 定义常数
m = 1;    % 质量
k = 1;    % 弹性系数
gamma = 0.1;  % 阻尼系数

% 定义初始值
x0 = 1;
y0 = 0;
vx0 = 0;
vy0 = 1;

% 定义时间间隔和离散点数
dt = 0.01;
t = 0:dt:100;

% 定义四阶龙格-库塔法求解函数
f = @(t,y) [y(3); y(4); -k/m*y(1) - gamma/m*y(3); -k/m*y(2) - gamma/m*y(4)];
[t, y] = ode45(f, t, [x0;y0;vx0;vy0]);

% 绘制位移和速度随时间变化的图像
subplot(2,1,1)
plot(t, y(:,1), t, y(:,2))
xlabel('t')
ylabel('x, y')
legend('x', 'y')

subplot(2,1,2)
plot(t, y(:,3), t, y(:,4))
xlabel('t')
ylabel('vx, vy')
legend('vx', 'vy')

这里，我们定义了质量 m、弹性系数 k、阻尼系数 gamma，以及初始位置和速度 x0、y0、vx0、vy0。然后，我们定义了时间间隔 dt 和离散点数 t，并使用 ode45 函数和定义好的微分方程 f 求解了离散解 y。

接下来，我们可以使用傅里叶变换函数 fft 计算离散解的功率谱：

% 计算功率谱
Fs = 1/dt;   % 采样率
L = length(y(:,1));  % 信号长度
f = Fs*(0:(L/2))/L;  % 频率向量

X = fft(y(:,1));
Pxx = abs(X/L).^2;
Pxx = Pxx(1:L/2+1);

Y = fft(y(:,2));
Pyy = abs(Y/L).^2;
Pyy = Pyy(1:L/2+1);

% 绘制功率谱图像
figure()
subplot(2,1,1)
plot(f, Pxx)
xlabel('f (Hz)')
ylabel('Pxx')

subplot(2,1,2)
plot(f, Pyy)
xlabel('f (Hz)')
ylabel('Pyy')

这里，我们首先计算了采样率 Fs、信号长度 L 和频率向量 f。然后，我们使用 fft 函数计算了 x(t) 和 y(t) 的傅里叶变换，并计算了功率谱 Pxx 和 Pyy。最后，我们使用 subplot 函数绘制了

该回答引用GPTᴼᴾᴱᴺᴬᴵ
根据你提供的代码，可能有以下几个问题：

1.傅里叶变换中的采样频率 f 和角频率 omega 的计算不一致，可能导致横坐标不正确。建议将计算采样频率的代码改为：

f = 1 / h; % 采样频率

2.傅里叶变换中的 delta 和 R1 的值未给出，无法确定 omega_0 的值。建议在代码中给出这些参数的值。

3.绘制功率谱时，应该将横坐标设置为频率，而非角频率。建议将绘图代码改为：

f = 1 / h; % 采样频率
f_Nyquist = f / 2; % Nyquist 频率
f_step = f / length(y(:,2)); % 频率步长
f_axis = 0:f_step:f_Nyquist; % 频率坐标轴
S_y = abs(fft(y(:,2))) .^ 2; % 功率谱
plot(f_axis, S_y(1:length(f_axis))); % 绘图

这样修改后，应该可以正确绘制功率谱了。


clear;
clc
fwind=0.01144;
R1=3.17;
R3=-1.61;
R5=-2.76;
g=9.81;
I=5.37*10^10;    
delta=1.47*10^8; 
r1=delta*g*R1/I;
r3=delta*g*R3/I;
r5=delta*g*R5/I;
 
%options=odeset('maxstep',0.01);%%时间步长调整
%[t,y]=ode45(@vdp1,[-30 30],[0.875;0],options);%%四阶龙哥库塔法求解微分方程组，y(:,1)即x的解，y(:,2)即y的解
odefun2 = @(t,y)[y(2);-(r1*y(1)+r3*y(1)^3+r5*y(1)^5)+fwind*(1-y(1)^2/2+y(1)^4/24)]; % 化二阶常微分方程为二阶常微分方程组
 
ts = -30; %计算开始时间
te = 30; %计算结束时间
h = 0.001;  %时间步长
t = ts:h:te; %计算数据点
y = zeros(length(t),2); 
y(1,:)=[0.875,0]; %参数初值
 
for n = 1:length(t)-1
k1 = odefun2(t(n),y(n,:)');
k2 = odefun2(t(n)+h/2,y(n,:)'+h/2*k1); %这里要特别注意同型向量相加
k3 = odefun2(t(n)+h/2,y(n,:)'+h/2*k2);
k4 = odefun2(t(n)+h,y(n,:)'+h*k3);
y(n+1,:) = y(n,:)+h/6*(k1+2*k2+2*k3+k4)';%四级四阶Runge-Kutta公式
end
%%傅里叶变换
N=length(y(:,2));
y_OMEGA=fftshift(fft(y(:,2)));
% f=(0:N-1)/(N*h);%%采样频率,这里不对，时间既然有负的，这里也得有，修改如下：
f = -1/(2*h):(1/(te-ts)):1/(2*h);
omega_0=(delta*g*R1/I)^0.5;  %问一个问题，这里这个是谱半径吧，没有2pi,f也不用乘
% omega=f*2*pi;
omega=f;
OMEGA=omega/omega_0; %求出相应横坐标Ω

%% 科普一下，频谱一般会加一个log(F)函数，或者10log10(F)
S_y=log((abs(y_OMEGA))); %Sy(Ω)

%% 右图处理，只画了正频率，因为是对称的，所以你不需要全部画，修改如下
OMEGA = OMEGA(ceil(N/2):end);
S_y = S_y(ceil(N/2):end);
plot(OMEGA,S_y,'b','Linewidth',2);
axis([0 1.5 0 inf]);legend('fwind=0.01144');

%% 要想更平滑，插值就行了，这用在发论文中可以，自己分析还是别了

运行结果如下

img