马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
function x=f(q)
b=0.1;
r=1.0;
a=0.3;
w=4:0.1:8;
Hlw=1/sqrt(8);
errorq=1;
q=1;
while errorq>0.00001
syms y;
b0=r*(2*a^2*int(1/sqrt(2*pi)*exp(-y^2/2),0,a/q)+2*a*q/sqrt(2*pi)*exp(-a^2/(2*q^2)));
b2=r*(4*a*int(1/sqrt(2*pi)*exp(-y^2/2),0,a/q)+4*q/sqrt(2*pi)*exp(-a^2/(2*q^2)));
b5=2*r*int(1/sqrt(2*pi)*exp(-y^2/2),0,a/q);
for j=1:41
for i=1:41
Hx1(i)=Hlw/(1-w(i)^2+sqrt(-1)*w(i)*(b+b2));
Hx11(i)=sqrt(-1)*w(i)*Hx1(i);
Hx2(i,j)=(b5*w(i)*w(j)*Hx1(i)*Hx1(j))/(1-(w(i)+w(j))^2+sqrt(-1)*(w(i)+w(j))*(b+b2));
Hx22(i,j)=sqrt(-1)*(w(i)+w(j))*Hx2(i,j);
end
end
Q=sum((abs(Hx11(i)))^2*0.1);
errorq=subs(abs(Q-q));
q=Q;
x=q;
end
上面是我编写的一个m函数,我想通过迭代求出q,(中间有好几步是对q没有关系的,但是我还要用到,所以没有删去)。运行之后给出的答案如下:ans =
4/5/(3969+(4/5+21560115664298739/11258999068426240*erf(24012833680751548415149939247787/30948500982134506872478105600*2^(1/2)+36993277947670772725057411478259/1980704062856608439838598758400*erf(3/20*2^(1/2))*pi^(1/2)+1394515762373819567048353910970363/1014120480182583521197362564300800*2^(1/2)*erf(3/20*2^(1/2))^2*pi)*2^(1/2)*pi^(1/2)+288230376151711744/28222125408961305/(3969+(571940004461227/43980465111040+21560115664298739/11258999068426240*erf(3/20*2^(1/2))*2^(1/2)*pi^(1/2))^2)*exp(-9/128*(3969+(571940004461227/43980465111040+21560115664298739/11258999068426240*erf(3/20*2^(1/2))*2^(1/2)*pi^(1/2))^2)^2))^2)。
我觉得肯定不对,所以请高手指点一下。不胜感激!(程序有点长,但是程序比较简单,所以希望高手不要嫌麻烦)
[ 本帖最后由 upc_uk 于 2007-9-5 08:49 编辑 ] |