请问这种涂改如何 做出来？

likefeng

大家好，我用ODE45解了一个微分方程组，现想用此方程组做出如附件的图，应该是幅频图。请问大家如何做？

likefeng

图中横坐标是W0/W1,
W0是方程中的W，方程组如下:

dx=[x(2);

0.002337*(F1-C12*x(4)-K11*x(1)-K12*x(3))-0.02*w1*x(2)-w1^(2)*x(1);

x(4);

0.00117357*(F2-C21*x(2)-K21*x(1)-K22*x(3))-0.02*w2*x(4)-w2^(2)*x(3)];

数值解程序如下：

function dx=gfz1(t,x,F1,F2,C12,C21,K11,K12,K21,K22,w1,w2)

s=1.6;b=1;w=4.8;w1=5.14;w2=9.49;

A=-0.3*w^2*cos(w*t).*cos(s*t)-s*(-0.3)*w*sin(s*t);

B=0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*cos(s*t);

C=-0.3*w^2*cos(s*t);

f1=48535.1+993.93*cos(s*t);

f2=-10426.4+(-177.19)*cos(s*t);

F1=f1+(-28510.61)*A.*b+(-6619.84)*B.*b^2+(-872.82)*C;

F2=f2+468.19*A*b+(-47391.25)*B*b^2+75.82*C;

C12=1078.62*(-0.3)*b*sin(s*t);

C21=1078.62*0.3*b*sin(s*t);

K1=b^2*0.09*w^(2)*sin(w*t).*sin(w*t);

K2=b*((-0.3)*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);

K3=b*(0.3*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);

K4=-b*(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*sin(s*t)-s^2;

K5=b^(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*cos(s*t).*cos(s*t);

K11=381.34*K1+81.22*K2+81.22*K3+17.72*K4+806.29*K5;

K12=-115.92*K1+517.92*K2+(-21.4)*K3+115.95*K4+(-206.29)*K5;

K21=-115.92*K1-21.42*K2+517.92*K3+115.95*K4+(-206.29)*K5;

K22=61.04*K1-114.93*K2-114.93*K3+784.94*K4+1373.31*K5;

dx=[x(2);

0.002337*(F1-C12*x(4)-K11*x(1)-K12*x(3))-0.02*w1*x(2)-w1^(2)*x(1);

x(4);

0.00117357*(F2-C21*x(2)-K21*x(1)-K22*x(3))-0.02*w2*x(4)-w2^(2)*x(3)];

T0=2*pi/w;t=0:T0/500:T0;

x0=[0;0;0;0];

[t,x]=ode45('gfz1',t,x0);

dujunmin

没有能太看懂上面的程序，感觉求出的应当是时域中的特征，如果是这样的话，那么首先要作傅立叶变换，才能得到幅频曲线。