matlab 方程不运算,MATLAB程序解方程组画图，运算好几天都没出来结果！ - 仿真模拟 - 小木虫 - 学术 科研 互动社区...

CODE:

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;