zhailiangjun 发表于 2009-8-19 16:14

求助 程序的简化

各位大侠,这是我的写的一个程序,主要的目的是把一个(N+1)^2*(N+1)^2方阵约化至(N+1)*(N+1)的方阵。但是里面有几层的循环,做起来很浪费时间,请大家帮忙看看有没有办法调整一下。

这个程序的思路其实是很简单的,rho是一个(N+1)^2*(N+1)^2的方阵,是由(N+1)^2的列矩阵和(N+1)^2行矩阵相乘得到的。rho2是一个(N+1)*(N+1)的方阵,是由rho约化得来的。a(程序中用a1和a2表示)和b分别是两组长度为(N+1)的基矢,具有相同的形式,每一组中的第n个矢量可以写成一个单位方矩阵的第n列元素。rho2中的元素,
             rho2(k1,k2)=\sum [ kron(a(k1),b(i)) ]' * rho *kron(a(k2),b(i))。
kron(a(k1),b(i)) 表示 a 这组基矢中的第 k1 个矢量和 b 中第 i 个矢量的直积,\sum是对b的所有基矢的求和。
就是这样了,请大家帮忙看看吧。多谢了。
function rho2=reduce(psi1,psi11,N)
format long
% global N
rho=psi1*psi11;
% psi1是一(N+1)^2维列矩阵,psi11为一(N+1)^2维行矩阵
rho2=zeros(N+1,N+1);
I=eye(N+1,N+1);
for k1=1:N+1
    a1=I(:,k1);
    for k2=1:N+1
      a2=I(:,k2);
      for k3=1:N+1
            b=I(:,k3);
            phi1=kron(a1,b);
            phi2=kron(a2,b);
            rho2(k1,k2)=phi1'*rho*phi2+rho2(k1,k2);
      end
    end
end
页: [1]
查看完整版本: 求助 程序的简化