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(-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 发表于 2007-9-5 09:35

原帖由 upc_uk 于 2007-9-5 08:47 发表 http://www.chinavib.com/forum/images/common/back.gif
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 发表于 2007-9-5 09:41

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

eight 发表于 2007-9-5 09:43

原帖由 xjzuo 于 2007-9-5 09:41 发表 http://www.chinavib.com/forum/images/common/back.gif
请将你的问题及公式用word上传一下,以便他人给出建议.
另:Hx1(i)*Hx1(j)处似乎不对,整个程序似乎也嫌累赘.

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

xjzuo 发表于 2007-9-5 09:48

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

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

xjzuo 发表于 2007-9-5 10:25

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

[ 本帖最后由 eight 于 2007-9-5 10:32 编辑 ]
页: [1]
查看完整版本: 能帮忙看看程序错在哪里么?