function dx=crackandnonlinerandnon(t,x,c1,w)
syms a k
e1=0.1;
R=0.008;
mu=0.1;
E=2*10^11;
L=0.5;
m=3.55;
%g=9.8;
n=6;
p=6;
q=6;
I=pi*R^4/4;
k0=48*E*I/L^3;
%st=m*g/k0;
%w0=sqrt(k0/m);
%c1=c/sqrt(m*k0);
%w1=w0*w;
r=(mu*(2-mu))^(1/2);
A1=R^2*(pi-acos(1-mu)+r*(1-mu));
e=2*R^3*r^3/(3*A1);
theta1=atan((e+R*(1-mu))/(R*r));
theta2=pi/2+acos(1-mu);
%感觉问题在下面这段
g1=(cos(a*theta2)-cos(a*theta1))*cos(a*t)/a^2;
g2=2*(0.8*theta2)^2*sin(pi-k*theta2*0.8)*sin(k*t)/(pi^2-k^2*theta2*0.8);
s1=symsum(g1,a,1,p);
s2=symsum(g2,k,1,q);
f1=(cos(t/2))^n;
f2=1/pi*((theta1+theta2)/2-2/(theta2-theta1)*s1);
%就是上面这段
Ix=pi*R^4/8-R^4/4*((1-mu)*(2*mu^2-4*mu+1)*r+asin(1-mu));
Iy=R^4/12*((1-mu)*(2*mu^2-4*mu-3)*r+3*asin(r));
Ixx=I-(Ix+A1*e^2)*f1;
Iyy=I+(Ix+A1*e^2)*f1-(Ix+Iy+A1*e^2)*f2;
Ixy=((Ix-Iy)/2+A1*e^2/2)*s2;
phi0=0;
dx=zeros(4,1);
dx(1)=x(2);
dx(2)=w^2+e1*cos(t+phi0)-2*c1*w*x(2)-w^2*(Ixx/I)*x(1)-w^2*(Ixy/I)*x(3);
dx(3)=x(4);
dx(4)=e1*sin(t+phi0)-2*c1*w*x(3)-w^2*(Ixy/I)*x(1)-w^2*(Iyy/I)*x(3);

下面是求解画图：
clear;clc;
tic;
c1=0.1;
j=1;
tend=linspace(0,2000,3500);
for w=0.1:0.05:5;
     x0=[0.01;0.01;0.01;0.01];
    options=odeset;options.RelTol=1e-4;
    [T,x]=ode45(@crackandnonlinerandnon,tend,x0,options,c1,w);
    ya=x(end-1500:end,:);
    R0=sqrt((ya(:,1).^2)+(ya(:,3).^2)); %求幅值
    [RM0,j]=max(R0);
     plot(w,RM0,'ro');
hold on;   
end

title('P-W图');
xlabel('回转速度W');
ylabel('振幅P');
time=toc;

3楼: Originally posted by 取经的僧人 at 2016-10-18 09:27:23
你花了多长时间仿真出来啊？我运行两三天都没有出来结果！跟电脑配置有关吗？我用的MATLAB2014！电脑64位8G内存！！！...