试过符号积分,没积出来。大家帮我想想办法吧。下面是我写的代码:
% ============================================================================
clear;clc
s0=0:0.05:0.3;
syms 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、SH、c
ch1=strcat(num2str(ss),'*cosh(dpdx*z/',num2str(ss),')/',num2str(v0));
sh1=strcat(num2str(ss),'*sinh(dpdx*z/',num2str(ss),')/',num2str(v0));
ch=subs(ch1,'dpdx',dpdx) %因为最终要对x积分,所以需要将dpdx关于x的表达式代入
sh=subs(sh1,'dpdx',dpdx) %因为最终要对x积分,所以需要将dpdx关于x的表达式代入
class(sh)
cH=inline(ch)
sH=inline(sh)
CH=quadl(cH,0,hc) %这一句有问题。这一句完成后,CH应该是x的函数。
SH=quadl(sH,0,hc)
c=(uc*CH-SH*(CH^2-SH^2+uc^2)^0.5)/(CH^2-SH^2)
%以下将CH、SH、c代入求被积函数fun
q=subs(strcat('dpdx*',num2str(hc),'+',num2str(ss),'*sinh(c)'),'dpdx',dpdx);
q=vpa(subs(q,'c',c));
fun=inline(q);
%以下积分求F2,进而求f
F1=quadl(fun,(-b),b)
F2=quadl(F1,(-L/2),L/2);
f(j)=double(F2/W);
end
ff(i,:)=f;
end
ff
plot(s0,ff,'-*')
% ================================================================================