求解2阶振动方程的matlab程序

star0304

u0=4*pi*1e-7;d0=0.5e-3;Fj=174.312;R=22.4e-3;L=56e-3;e=aqrt(x^2+y^2);
E=e/d0;w=60000*2*3.14/60; t=0:0.01：10;r=w*t;
A0=u0/d0/sqrt(1-E^2);
A1=2*u0/d0/sqrt(1-E^2)*(1-sqrt(1-E^2)/E);
A2=2*u0/d0/sqrt(1-E^2)*(1-sqrt(1-E^2)/E)^2;
A3=2*u0/d0/sqrt(1-E^2)*(1-sqrt(1-E^2)/E)^3;
f1=R*L*pi/4/u0*Fj^2*(2*A0*A1+A1*A2+A2*A3);
f2=R*L*pi/4/u0*Fj^2*(A0*A1+1/2*A1*A2+1/2*A2*A3);
f3=R*L*pi/4/u0*Fj^2*(A0*A3+1/2*A1*A2);
Fx=f1*cos(r)+f2*cos(2*w*t-r)+f3*cos(2*w*t-3*r);
Fy=f1*sin(r)+f2*sin(2*w*t-r)-f3*sin(2*w*t-3*r);
%质量、刚度矩阵是用有限元算的90x90的矩阵（程序以打包发送）；阻尼矩阵为C=0.05*M+0.02*K;
a Ω为常数 要求解这个方程​ 的matlab程序

追风少年-阿杰

我现在也要编程求解两自由度弹簧振子系统的振幅随激励频率的变化及相关参数的影响关系。请问你编出来了吗？