能帮忙看看程序错在哪里么？

upc_uk

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 编辑 ]

eight

原帖由 upc_uk 于 2007-9-5 08:47 发表

                               
登录/注册后可看大图

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(- ...

请先看看各个置顶帖，然后自己认真考虑一下你的帖子是否表述清楚了

xjzuo

请将你的问题及公式用word上传一下，以便他人给出建议．
另：Hx1(i)*Hx1(j)处似乎不对，整个程序似乎也嫌累赘．

eight

原帖由 xjzuo 于 2007-9-5 09:41 发表

                               
登录/注册后可看大图

请将你的问题及公式用word上传一下，以便他人给出建议．
另：Hx1(i)*Hx1(j)处似乎不对，整个程序似乎也嫌累赘．

呵呵，还是 xjzuo 版主你认真，我一看到这类帖子就不会细看了（光有代码，没有说明，至少这个代码是干什么的都没有说）

xjzuo

没办法,有时候作为版主, 总是不得不重复提醒版友"应该首先将问题讲清楚,别人才能给出好建议".

其实这只是一个简单的迭代求解问题, 由于他没有贴公式, 不好检查其错误之处,更不用说给出简洁的求解办法了------有时候猜测总是很累的.

xjzuo

我修改程序计算了一下, 结果为1.9998e-004, 不知道是否是LZ所求？

[ 本帖最后由 eight 于 2007-9-5 10:32 编辑 ]