|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
function w=kx2(x)
global D1 D2 a b k1 k2 tao1 tao2 lamda p q h g1 A B G G1 xs H;
w=H/(D2*(G1^2)+2*lamda*sqrt(D1*D2)*G1/B+D1);
global D1 D2 a b k1 k2 tao1 tao2 lamda p q h g1 A B G G1 xs H;
D1=0.1;
D2=0.2;
a=5.6;
b=0.4;
k1=10;
k2=1.0;
tao1=0.5;
tao2=0.25;
lamda=0.4;
dx=0.16;
fid=fopen('99.dat','wt')
for ii=1:25
x=0.01+dx*ii;
xx(ii)=x;
p=k2-(a+b)^2/(3*k1^2);
q=-2*(a+b)^3/(27*k1^3)+k2*(a+b)/(3*k1)-b;
A=1-tao1.*(2*a*k1.*x./((x.^2+k1).^2)-k2)+tao1*2*k1*(a.*x.^2./(x.^2+k1)-k2.*x+b)./(x.*(x.^2+k1));
h=a.*(x.^2)./(x.^2+k1)-k2.*x+b;
g1=x.^2./(x.^2+k1);
xs=(-q/2+sqrt((q/2)^2+(p/3)^3))^(1/3)+(-q/2-sqrt((q/2)^2+(p/3)^3))^(1/3);
B=1-tao2*(2*a*k1*xs/((k1+xs)^2)-k2);
H=h./A;
G1=g1./A;
G=sqrt(D2*(G1^2)+2*lamda*sqrt(D1*D2)*G1/B+D1);
m=quad(@kx2,0,x);
yy(ii)=exp(m-log(abs(G)));
end
sump=sum(yy);
for jj=1:25
yy(jj)=yy(jj)/sump;
Q=yy(jj);
fprintf(fid,'%f\n',Q/dx);
end
plot(xx,yy/dx);
fclose(fid);
hold on;
出错提示:
??? Attempted to access y(7); index out of bounds because numel(y)=1.
Error in ==> quad at 71
if ~isfinite(y(7))
Error in ==> xm2 at 27
m=quad(@kx2,0,x);
新手上路请多指教,谢谢!
|
|