矩阵计算的问题
有一列矩阵，它的第k个元素是这样写的，A(k,1)=^(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)=^(1/2)*(Z')^(N-k)*(1-Z*Z')^(1/2*k);
这样的两个矩阵相乘后B*A可以得到，
B*A=\sum *(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 的帖子
嗯，谢谢你的回复，你说的这个，我是特意考虑了的，取了共轭了，所以在N和P的值较小的时候是等于1的，但是大了就不行了。 阶乘函数的值比较大，大数+-小数容易吃掉小数，这是数值计算与生俱来的问题。 可以适当考虑将大的因子提出来，作点变形，保证两个相加减数的量级接近 试比较，理论结果真是1吗？% clc
% clear
function shiyan
p=10;
q=1;
k=8;
N=10;
Z=(q+i*p)/(2*N)^(1/2);
=myA(k,Z,N);
disp(B*A)
disp(sum(C))
function =myA(x1,x2,x3)
y1=zeros(x1,1);
y2=zeros(1,x1);
y3=zeros(1,x1);
for j=1:x1
y1(j)=sqrt(nchoosek(x3,j))*x2^(x3-j)*(1-x2'*x2)^(1/2*j);
y2(j)=sqrt(nchoosek(x3,j))*(x2')^(x3-j)*(1-x2'*x2)^(1/2*j);
y3(j)=nchoosek(x3,j)*(x2'*x2)^(x3-j)*(1-x2'*x2)^(j);
end结果：
2.8297e+006 +9.5367e-007i
2.8297e+006
[ 本帖最后由 friendchj 于 2009-7-18 22:01 编辑 ]
回复 5楼 friendchj 的帖子
嗯，这个公式(a+b)^n=\sum *a^(n-k)*b(k)
是不错的，当a+b=z'z+(1-z'z)=1的时候，那么这个结果就应该是等于1的啊。
谢谢你的计算，nchoosek这个函数非常好用。谢谢你。
回复 地板 VibrationMaster 的帖子
非常感谢你的建议谢谢你的留言，利用了5楼 friendchj的程序计算的时候，可以发现，当p的值取得很小的时候，例如令p=0，这个时候将N的值取得很大，计算出来的结果也是等于1的。但是当p的值取得稍微大些的时候，例如大于10的时候，结果就会有很大的差别了。阶乘函数的结果是和N的值相联系在一起的，所以在这个计算中，可能N的大小并不是一个影响很大的因素，重要还是在与z=(q+i*p)/(2N)^(1/2)这个地方，不知道怎么处理才好。
[ 本帖最后由 ChaChing 于 2010-4-25 16:34 编辑 ] 看了你的帖子，仔细分析了一下，如果要使结果为1，则n和k要满足k=n+1,例如(a+b)^2有3(k)项。
如果k小于n+1，则结果不是所有项的和，也不可能等于1了。
function shiyan
p=10;
q=1;
k=11;
N=10;
Z=(q+i*p)/(2*N)^(1/2);
=myA(k,Z,N);
disp(B*A)
disp(sum(C))
function =myA(x1,x2,x3)
y1=zeros(x1,1);
y2=zeros(1,x1);
y3=zeros(1,x1);
for j=1:x1
y1(j)=sqrt(nchoosek(x3,j-1))*x2^(x3-j+1)*(1-x2'*x2)^(1/2*(j-1));
y2(j)=sqrt(nchoosek(x3,j-1))*(x2')^(x3-j+1)*(1-x2'*x2)^(1/2*(j-1));
y3(j)=nchoosek(x3,j-1)*(x2'*x2)^(x3-j+1)*(1-x2'*x2)^(j-1);
end
结果：
1.0000 + 0.0000i
1.0000
回复 9楼 friendchj 的帖子
嗯，谢谢你那么热心的建议，受益匪浅。使用你这个程序计算的话还是只能在p<=10的这个范围之内很有效果，就是等于1，但是我要计算的东西一般是在P=200附近的，这样的话还是有很大的误差的。谢谢你。回复 10楼 zhailiangjun 的帖子
P很大时，需要N较小时结果才能正确，例如：p=200，N=3时结果为：0.9984 + 0.0002i
1.0001
如果P很大，N也很大，计算中会出现(z)^N这种形式，可能会导致结果不正确。 1-Z'*Z 当P,N很大时会吃掉1 。必须改造原表达式
回复 11楼 friendchj 的帖子
呵呵，是的啊，但是我要计算时N的取值大概要在40左右的，所以这个就会误差很大了。呵呵，所以还是要想办法改写程序。谢谢你参与讨论。:lol回复 12楼 VibrationMaster 的帖子
是的，确实问题在这个1-Z'*Z 表达式上。谢谢你的讨论。:loveliness:回复 9楼 friendchj 的帖子
你好，关于这个小的矩阵的运算的时候，当 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
这样改写之后，会在一定程度上校正了矩阵，但是仍然不能得到很正确的结果。
请大家再帮忙看看，多谢了。
页:
[1]