用MATLAB编写s函数举例,S函数编写 - 仿真模拟 - MATLAB - 小木虫论坛-学术科研互动平台...

如何把下列方程组程序编辑到S函数？

function dx=crackandnonlinerandnon(t,x,c1,w)

e1=0.1;

R=0.008;

mu=0.1;

E=2*10^11;

L=0.5;

n=10;

p=10;

q=10;

I=pi*R^4/4;

K0=48*E*I/L^3;

r=sqrt(mu*(2-mu));

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);

a=1:p;

g1=(cos(a*theta2)-cos(a*theta1)).*cos(a*w*t)./a.^2;

s1=cumsum(g1);

ss1=s1(length(t));

k=1:q;

g2=2*0.8*theta2*sin(pi-k*theta2*0.8).*sin(k*w*t)./(pi^2-k.^2*(theta2*0.8)^2);

s2=cumsum(g2);

ss2=s2(length(t));

f1=(cos(w*t/2))^n;

f2=-1/pi*((theta1+theta2)/2-2*ss1/(theta2-theta1));

Ix=pi*R^4/8-R^4*((1-mu)*(2*mu^2-4*mu+1)*r+asin(1-mu))/4;

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)*ss2;

Kxx=48*E*Ixx/L^3;

Kyy=48*E*Iyy/L^3;

Kxy=48*E*Ixy/L^3;

phi0=0;

ep1c=0.1;%非线性参量

ep1s=0.1;

ep3c=0;

ep3s=0;

ep30=ep1c+ep3c;

ep21=ep1s+3*ep3s;

ep12=ep1c-3*ep3c;

ep03=ep1s-ep3s;

be0=0.1;

be2c=0;

be4c=0;

be2s=0;

be4s=0;

be40=be0+be2c+be4c;

be31=2*be2s+4*be4s;

be22=2*be0-6*be4c;

be13=2*be2s-4*be4s;

be04=be0-be2c+be4c;

Nx=3*ep30*x(1)^2+2*ep21*x(1)*x(3)+ep12*x(3)^2+4*be40*x(1)^3 ...

+3*be31*x(1)^2*x(3)+2*be22*x(1)*x(3)^2+be13*x(3)^3;

disp('进入Ny');

Ny=ep21*x(1)^2+2*ep12*x(1)*x(3)+3*ep03*x(3)^2+be31*x(1)^3 ...

+2*be22*x(1)^2*x(3)+3*be13*x(1)*x(3)^2+4*be04*x(3)^3;

dx=zeros(4,1);

dx(1)=x(2);

dx(2)=1+e1*w^2*cos(w*t+phi0)-c1*x(2)-(Kxx/K0)*x(1)-(Kxy/K0)*x(3)-Nx;

dx(3)=x(4);

dx(4)=e1*w^2*sin(w*t+phi0)-c1*x(4)-(Kxy/K0)*x(1)-(Kyy/K0)*x(3)-Ny;