snowfish 发表于 2008-11-11 15:45

大家能看看问题出在哪吗

下面是我写的一个完整的小程序，大家能不能帮我试下，为什么delta循环到0.21以后就没有值了呢，这个i2的值应该是在i3的附近，比i3大一点点的，在fzero中，我取的初始值就为I3，前面几个数据都挺好的，到了0.21以后i2就是nan了？很奇怪，我不知道问题出在哪啊？大家能不能运行看下， 非常感谢哦~:loveliness:
format long
tic;
k=0.01;
ome=2.5;
N=1000;
g20=0.1;
epson=0.1;
gam1=0.5*4;
gam2=0.5;
count_1=0;
delta=;
for xx=delta;
   normaldelta=xx
cos2=1/2+xx/(2*sqrt(1+xx^2));
sin2=1-cos2;
A1=ome*sqrt(cos2);
I1=A1*A1;

G2=g20*sqrt(cos2-sin2);
cosh=g20*cos2/G2;
sinh=g20*sin2/G2;
sin4=sin2*sin2;
cos4=cos2*cos2;
emr=cosh-sinh;
em2r=emr^2;

M=cosh*sinh;
N0=sinh*sinh;

E1=0;
R1=gam1*cos2;
E2=gam2*sin4;
R2=gam2*cos4;
E3=gam1*sin2;
R3=0;
gamph=gam2*cos2*sin2;
      R12=R1+E1+R3+R2+gamph;
      R23=R1+R2+E2+E3+4*gamph;
      R13=E1+R3+E2+E3+gamph;

      U0=-(E1+R3)*(R2+E2+E3)-R1*(E2+E3+R3)+E3*(R3-R2);
      U1=R2+2*E2+2*E3+R3;
      U2=E2+E3+2*R3-R2;
      U3=2*E1+2*R3+R1+E3;
      U4=E1+R3-R1+2*E3;
      U5=3*R13;
      A11=-U1;
      A22=-U3;
      A12=U2+U4-U5;
      AA1=U0*R12-U1*R13*R23;
      AA2=U0*R23-U3*R13*R12;
      A0=U0*R12*R23*R13;
      B2=A12*I1+AA2;
      C2=A11*I1*I1+AA1*I1+A0;
      V0=-(E2+E3)*(R1-E1-R3)-E3*(R2+R3);
      V1=-(R1+E3-E1-R3);
      V2=E2+E3-R2-R3;
      W0=R3*(R1+E3)+(E1+R3)*(R2-E2-E3);
      B1=G2*G2*N*W0;
      C1=G2*G2*N*((V0-V2*R13)*I1+W0*R13*R12);
      A=A22;
      B=B2-B1/k;
      C=C2-C1/k;

F=B*B-4*A*C;
if F>0
   I0=(-B-sqrt(F))/2/A;
else
    I0=0;
end
   ifI0>0
   I3=I0
else
   I3=0.000000001
end

ep=epson*sqrt(cos2-sin2) ;
u=inline('ep/sqrt(y)+(B1*y+C1)/(A*y^2+B2*y+C2)-k','y','ep','B1','C1','A','B2','C2','k');
I2=fzero(u,I3,[],ep,B1,C1,A,B2,C2,k)

Y=I2;
count_1=count_1+1;
Y0(count_1)=Y;
end
Y0;

delta=delta(~isnan(Y0));
Y0(isnan(Y0))=[];
data=
save('f:\coherent123.dat','data','-ASCII')
plot(delta,Y0)

ch_j1985 发表于 2008-11-11 20:56

回复 楼主 snowfish 的帖子

请问LZ这是关于哪方面知识的程序？
头都大了，没有头绪……

ChaChing 发表于 2008-11-11 23:04

下午有试着看下, 真的很难有耐心
建议楼主稍微解释下, 别人较有可能帮忙!

ChaChing 发表于 2008-11-12 11:30

回复 楼主 snowfish 的帖子

唉! 真是苦功! 我试run下并进入debug模式协助侦错!
发现楼主inline给定的函数并非for all x, 也就是负值的x函数值会变复数(complex)!
我试着改fzero成 I2=fzero(u,,[],ep,B1,C1,A,B2,C2,k)
可以跑出东东, 但我不知道是否是楼主要的!

snowfish 发表于 2008-11-12 15:40

谢谢大家了，我也发现了，是fzero函数里面的初始搜索值的问题，那个值不晓得怎么给好，给的不合适的话，就只有某一段范围内有值。
大家有没有什么好的方法啊？

snowfish 发表于 2008-11-12 15:46

比如，我画 delta=范围之内的值，而方程的根I2大概分布在0-200之间，如果我给I2=fzero(u,100,[],ep,B1,C1,A,B2,C2,k) 初始搜索值为100的话，前面一小段，delta=范围内就没有12的值了。。:@L

ChaChing 发表于 2008-11-12 16:27

个人以为不同函数求零解, 本应有其不同的工程应用范围, 很难一体适用的
可以根据工程经验或绘出函数图形, 帮助给定合适范围, 避免数值求解方法发散
页: [1]
查看完整版本: 大家能看看问题出在哪吗