矩阵计算的问题

zhailiangjun

有一列矩阵，它的第k个元素是这样写的，
A(k,1)=[N!/(k!*(N-k)!)]^(1/2)*(Z)^(N-k)*(1-Z'*Z)^(1/2*k);
其中，N为参数，Z是一个复数，Z=(q+i*p)/(2*N)^(1/2)，其中，q，p均为实数。N!表示N的阶乘，k!,(N-k)!也是表示阶乘，Z'表示复数Z的共轭。
有一个行矩阵，它的第k个矩阵元写为，
B(1,k)=[N!/(k!*(N-k)!)]^(1/2)*(Z')^(N-k)*(1-Z*Z')^(1/2*k);
这样的两个矩阵相乘后B*A可以得到，
B*A=\sum [N!/(k!*(N-k)!)]*(Z'*Z)^(N-k)*(1-Z'*Z)^(k)
      =(Z'*Z+1-Z'*Z)^(N)
      =1
这个是理论结算的结果是这样，B*A=1。但是在用matlab编程计算时出现了这样一个情况，就是当我选择的p的值以及N的值较大的时候，B*A就不再等于1了，而是一个相当大的值，这个情况该怎么办呐。请大家多多帮忙。感激不尽。

VibrationMaster

取转置的时候看看是否取共扼了

zhailiangjun

嗯，谢谢你的回复，你说的这个，我是特意考虑了的，取了共轭了，所以在N和P的值较小的时候是等于1的，但是大了就不行了。

VibrationMaster

阶乘函数的值比较大，大数+-小数容易吃掉小数，这是数值计算与生俱来的问题。 可以适当考虑将大的因子提出来，作点变形，保证两个相加减数的量级接近

friendchj

试比较，理论结果真是1吗？

1. % clc
   
2. % clear
   
3. function shiyan
   
4. p=10;
   
5. q=1;
   
6. k=8;
   
7. N=10;
   
8. Z=(q+i*p)/(2*N)^(1/2);
   
9. [A,B,C]=myA(k,Z,N);
   
10. disp(B*A)
    
11. disp(sum(C))
    
12. 
    
13. function [y1,y2,y3]=myA(x1,x2,x3)
    
14. y1=zeros(x1,1);
    
15. y2=zeros(1,x1);
    
16. y3=zeros(1,x1);
    
17. for j=1:x1
    
18. y1(j)=sqrt(nchoosek(x3,j))*x2^(x3-j)*(1-x2'*x2)^(1/2*j);
    
19. y2(j)=sqrt(nchoosek(x3,j))*(x2')^(x3-j)*(1-x2'*x2)^(1/2*j);
    
20. y3(j)=nchoosek(x3,j)*(x2'*x2)^(x3-j)*(1-x2'*x2)^(j);
    
21. end

复制代码

结果：
  2.8297e+006 +9.5367e-007i
  2.8297e+006

[ 本帖最后由 friendchj 于 2009-7-18 22:01 编辑 ]

zhailiangjun

嗯，这个公式
(a+b)^n=\sum [n!/((k!)*(n-k)!)]*a^(n-k)*b(k)
是不错的，当a+b=z'z+(1-z'z)=1的时候，那么这个结果就应该是等于1的啊。
谢谢你的计算，nchoosek这个函数非常好用。谢谢你。

zhailiangjun

非常感谢你的建议
谢谢你的留言，利用了5楼 friendchj的程序计算的时候，可以发现，当p的值取得很小的时候，例如令p=0，这个时候将N的值取得很大，计算出来的结果也是等于1的。但是当p的值取得稍微大些的时候，例如大于10的时候，结果就会有很大的差别了。阶乘函数的结果是和N的值相联系在一起的，所以在这个计算中，可能N的大小并不是一个影响很大的因素，重要还是在与z=(q+i*p)/(2N)^(1/2)这个地方，不知道怎么处理才好。

[ 本帖最后由 ChaChing 于 2010-4-25 16:34 编辑 ]

friendchj

看了你的帖子，仔细分析了一下，如果要使结果为1，则n和k要满足k=n+1,例如(a+b)^2有3(k)项。
如果k小于n+1，则结果不是所有项的和，也不可能等于1了。

1. function shiyan
   
2. p=10;
   
3. q=1;
   
4. k=11;
   
5. N=10;
   
6. Z=(q+i*p)/(2*N)^(1/2);
   
7. [A,B,C]=myA(k,Z,N);
   
8. disp(B*A)
   
9. disp(sum(C))
   
10. function [y1,y2,y3]=myA(x1,x2,x3)
    
11. y1=zeros(x1,1);
    
12. y2=zeros(1,x1);
    
13. y3=zeros(1,x1);
    
14. for j=1:x1
    
15. y1(j)=sqrt(nchoosek(x3,j-1))*x2^(x3-j+1)*(1-x2'*x2)^(1/2*(j-1));
    
16. y2(j)=sqrt(nchoosek(x3,j-1))*(x2')^(x3-j+1)*(1-x2'*x2)^(1/2*(j-1));
    
17. y3(j)=nchoosek(x3,j-1)*(x2'*x2)^(x3-j+1)*(1-x2'*x2)^(j-1);
    
18. end

复制代码

结果：
   1.0000 + 0.0000i
    1.0000

zhailiangjun

嗯，谢谢你那么热心的建议，受益匪浅。使用你这个程序计算的话还是只能在p<=10的这个范围之内很有效果，就是等于1，但是我要计算的东西一般是在P=200附近的，这样的话还是有很大的误差的。谢谢你。

friendchj

P很大时，需要N较小时结果才能正确，例如：p=200，N=3时结果为：
0.9984 + 0.0002i

    1.0001
如果P很大，N也很大，计算中会出现(z)^N这种形式，可能会导致结果不正确。

VibrationMaster

1-Z'*Z 当P,N很大时会吃掉1 。必须改造原表达式

zhailiangjun

呵呵，是的啊，但是我要计算时N的取值大概要在40左右的，所以这个就会误差很大了。呵呵，所以还是要想办法改写程序。谢谢你参与讨论。:lol

zhailiangjun

是的，确实问题在这个1-Z'*Z 表达式上。谢谢你的讨论。:loveliness:

zhailiangjun

你好，关于这个小的矩阵的运算的时候，当 1-z'*z<0时，就是说q^2+p^2>2N时就会出现较大的误差，一个原因就是程序中的 (1-x2'*x2)^(1/2*(j-1)) 这一项，当1-x2'*x2<0时，这一项应该是一个纯虚数，例如(-50)^(3/2)应该等于 i*3.535533905932738e+002，但是在matlab的运算中得到的结果是这样的，(-50)^(3/2)
ans =
    -6.494670421766197e-014 -3.535533905932737e+002i
不仅有个很小的实部存在，而且虚部的符号也是相反的，会造成很大的误差。可以改写成
sqrt((1-z'*z)^(j)).例如：
sqrt((-50)^3)
ans =
                  0 +3.535533905932738e+002i
这样改写之后，会在一定程度上校正了矩阵，但是仍然不能得到很正确的结果。
请大家再帮忙看看，多谢了。