今天看了您的示例帖受益匪浅。后来,又把自己的程序改了一下,好了一些。不过还是有点问题。刚学不久,很多问题都不懂。请版主多多指教。
clear;clc
s0=0:0.05:0.3;
syms l1 x z dpdx
n2=100; l=16.52; L=l*0.001; R=0.015; E=2.2637e11;
v0=0.03527432; vi=2.19e-8; ss=3e6;
ff=zeros(6,length(s0));
for i=1:6
W=i*1000;
f=[];
for j=1:length(s0)
s1=s0(j);
u2=n2*pi*R/30; u1=((2-s1).*u2)./(2+s1); u=(u1+u2)./2; uc=u2-u1;
hmin=(2.65*vi^0.54*(v0*u).^0.7*R^0.43*L^0.13/E^0.03)./W.^0.13;
hc=hmin*4/3;
b=(8*W*R/pi/L/E).^(0.5);
ph=2*W/L/b/pi;
p=ph*(1-x.^2/b^2).^(1/2);
dpdx=vpa(diff(p));
ch=strcat(num2str(ss),'*cosh(dpdx*z/',num2str(ss),')/',num2str(v0));
sh=strcat(num2str(ss),'*sinh(dpdx*z/',num2str(ss),')/',num2str(v0));
cH=vpa(subs(ch,'dpdx',dpdx));
sH=vpa(subs(sh,'dpdx',dpdx));
CH=int(cH,z,0,hc);
SH=int(sH,z,0,hc);
c=vpa((uc*CH-SH*(CH^2-SH^2+uc^2)^0.5)/(CH^2-SH^2));
q=subs(strcat('dpdx*',num2str(hc),'+',num2str(ss),'*sinh(c)'),'dpdx',dpdx);
q=vpa(subs(q,'c',c));
fun=inline(q);
F1=quadl(fun,(-b),b)
F2=int(F1,l1,(-L/2),L/2);
f(j)=double(F2/W);
end
ff(i,:)=f;
end
ff
plot(s0,ff,'-*') |