线性方程组的求解结果不一样的问题

回忆的路上

我编了个程序想验证下线性方程组求解
但是结果不对
需要解的方程组如下所示
要求解P1,P2,P3,P4
H(w1),H(w2),H(w3),H(w4)
是由要求解的已知的公式求出来的
应该解方程结果和我所用的P1,P2,P3,P4是一样的
但是结果不一样
程序如下:
clc;
clear;
N0 = 2048;
fs = 10e6;
f0 = fs/N0;
Omega0 = 2*pi*f0;
Omega  = [0:N0-1]*Omega0;
P    = [0.575  0.362  0.057  0.006];
tau  = [    0 0.2e-6 0.4e-6 0.6e-6];
H_calcu = P(1)*exp(-j*Omega*tau(1)) ...      % 我用这个公式生成的数据
        + P(2)*exp(-j*Omega*tau(2)) ...
        + P(3)*exp(-j*Omega*tau(3)) ...
        + P(4)*exp(-j*Omega*tau(4));
H_calcu_abs = abs(H_calcu);
seq = 1:4;                                                  % 因为要求解P(1)...P(4)四个参数,所以选择4个数据
Omega_x = Omega(seq);                             % 对应公式中用到的Omega也取对应的4个数据
H       = H_calcu(seq);
% 下面解方程组,计算结果与我所使用的P(1)...P(4)不同
W_matrix = [exp(-j*Omega_x(1)*tau(1)) exp(-j*Omega_x(1)*tau(2)) exp(-j*Omega_x(1)*tau(3)) exp(-j*Omega_x(1)*tau(4)); ...
            exp(-j*Omega_x(2)*tau(1)) exp(-j*Omega_x(2)*tau(2)) exp(-j*Omega_x(2)*tau(3)) exp(-j*Omega_x(2)*tau(4)); ...
            exp(-j*Omega_x(3)*tau(1)) exp(-j*Omega_x(3)*tau(2)) exp(-j*Omega_x(3)*tau(3)) exp(-j*Omega_x(3)*tau(4)); ...
            exp(-j*Omega_x(4)*tau(1)) exp(-j*Omega_x(4)*tau(2)) exp(-j*Omega_x(4)*tau(3)) exp(-j*Omega_x(4)*tau(4)) ];
H_tau    = H;
Pl = (inv(W_matrix)*H_tau')';
不知道这是怎么回事?高手帮帮忙吧
谢谢

[ 本帖最后由 eight 于 2007-9-12 11:21 编辑 ]

回忆的路上

这个是方程

xjzuo

问题讲得较清楚, 值得鼓励.
－－－尽量不要用inv, 而用斜杠求解．
提示: 分清楚左除和右除，你的问题就解决了．

回忆的路上

谢谢xjzuo
我用斜杠也求过,结果是一样的
现在问题可能在于这个线性方程组的系数矩阵不是满秩的
或者本来是满秩的,matlab的精度不够
我用rank(matrix)求出来有时候(选择不同的w值)是1,2,3
而不是4,我怎么样才能选择一个满秩的系数矩阵?

回忆的路上

起初以为是矩阵奇异的问题
但是刚刚看了些资料
那个系数矩阵是个vandermonde
而且也满足vandermonde矩阵非奇异的条件啊
这是怎么回事?

回忆的路上

刚试了下W_matrix这个矩阵肯定是满秩的
大家帮我看看问题在哪吧?
而且取不同的四个Omega值做计算
得到的Pl结果还不一样
今天都第二天了,我不想被这个小问题挡脚啊
谢谢大家

[ 本帖最后由 回忆的路上 于 2007-9-13 11:33 编辑 ]

xjzuo

按我的提示去做, 相信就能自己找到原因.
另: 我已经算过了, "分清楚左除和右除"后, 没有任何问题.

回忆的路上

非常谢谢您耐心的解答
我在matlab的help中复制了一段话:
If A is a square matrix, A\B is roughly the same as inv(A)*B, except it is computed in a different way. If A is an n-by-n matrix and B is a column vector with n elements, or a matrix with several such columns, then X = A\B is the solution to the equation AX = B computed by Gaussian elimination with partial pivoting (see Algorithm for details). A warning message is displayed if A is badly scaled or nearly singular.
我的问题应该就是这个问题:
W_matrix * Pl = H_tau;     (1)式
其中W_matrix是一个4*4的矩阵,Pl是一个4*1的矩阵,H_tau也是一个4*1的矩阵
我改过程序了,求解Pl用下式:
Pl = (W_matrix\H_tau')';    (2)式
由(2)式计算得到的结果代入(1)式,结果只有1e-15量级的误差,应该是matlab计算误差
这说明(1)式的计算是对的
但是计算得到的Pl并不等于P啊?
我得到的Pl=[2.7127 - 0.027131i       -3.5808 + 0.081386i       2.5175 - 0.08138i      -0.64941 + 0.027124i];
而我用来计算的P = [0.575  0.362  0.057  0.006];
能指教下我哪个细节做错了吗?
谢谢!

[ 本帖最后由 回忆的路上 于 2007-9-13 15:40 编辑 ]

xjzuo

唉! 我已经强调两遍了: "分清楚左除和右除", 也不知道你有没有思考一下我的提示.
------你自己分别试试左除和右除不就清楚了? 也用不着给我发短消息...

回忆的路上

问题解决了,谢谢你了,呵呵
问题不在于左除右除
而是复数矩阵行向量转列向量应该用H.'    而不是H'
后面提到的那种是共轭转置

xjzuo

看来你很坚持,不过也不一定是坏事.
但有时候能接受别人的建议也许会使你少走一些弯路:
(虽然你上面的做法数值上是等价的)
%%%%%%%%%%%%%%%%%
P1=H_tau/W_matrix

P1 =
   0.5750 - 0.0000i   0.3620 + 0.0000i   0.0570 - 0.0000i   0.0060 + 0.0000i

回忆的路上

谢谢您的解答
另外，我想就左除右除这个问题再跟您讨论一下
我先简要述说一下两种运算（X为待求量）：
如果方程为：AX=B，那么如果有解的话，解为：X=A\B，这就是左除
如果方程为：XA=B，那么如果有解的话，解为：X=B/A，这就是右除
具体解释可以在matlab中mldivide和mrdivide的help中找到
我遇到的问题应该用左除来解，即：X=A\B
下面我来解释为什么你用右除的算式——P1=H_tau/W_matrix，也可以得到数值上正确的答案
首先定义转置算子为T，不论该算子为共轭转置还是非共轭转置，只要矩阵的维数和秩满足要求即可，
那么下面的推论成立：
如果A*B=C，那么T(A*B)=T(B)*T(A)=T(C)，其中A，B，C都为能满足该等式运算维数和秩的矩阵
如果考虑方程组：AX=B，那么根据上面推论可以得到：T(A*X)=T(X)*T(A)=T(B)
这个时候用矩阵右除可以得到解：T(X)=T(B)/T(A)
你所用的矩阵右除的计算式：P1=H_tau/W_matrix中，
Pl就相当于T(X)，H_tau就相当于T(B)，W_matrix就相当于T(A)
在我遇到的问题中，T(X)和T(B)都是一维向量（矩阵），所以对他们取T算子运算都只是行/列向量的转变，或者共轭转置，而对于T(A)，在你的矩阵右除计算式中应该是T(W_matrix)，即应该是矩阵W_matrix的转置（或者共轭转置），我认为你式子的问题就在这里，数值计算结果正确，并不是因为你的式子正确，而是因为在解这个等式的过程中，不管取W_matrix，还会T(W_matrix)，都会得到近似相等的结果，这个结论你可以改下程序来验证。

xjzuo

你也不傻啊-----你讨论一下行向量的情形呢?
列向量的情形如你所言----这也是我为什么说"看来你很坚持,不过也不一定是坏事".
-----看来我不应该参与这个讨论,好自为之吧...