| 查看: 1227 | 回复: 0 | ||
[求助]
S函数编写
|
|
如何把下列方程组程序编辑到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; |
回复此楼» 猜你喜欢
求各位大神看下
已经有27人回复
-
影响面上的因素
已经有10人回复
-
咨询面上基金
已经有5人回复
-
8月时间戳变的，举个手。玩一下，释放压力
已经有11人回复
-
一个有机合成实验室都需要哪些设备？
已经有11人回复
-
面上再次挂了，太难了，躺也躺不了，倦也卷不过，小学校之殇！
已经有30人回复
-
2027年申博
已经有4人回复
-
好消息？这个有何含义？？？
已经有10人回复
-
【2027博士申请】纳米药物递送方向
已经有4人回复
-
UV压敏胶开发
已经有7人回复