matlab中十阶单位矩阵,用matlab或其他软件编程求解几十阶以上的矩阵

引用回帖:

virtualzx at 2012-11-02 11:24:59

我有点迷惑，你说让这个矩阵等于0，什么意思？是指矩阵的行列式等于0吗？那样的话不是一元一次方程，100个原子就是一个100次方程

是这样子的，我把程序贴出来syms m k mu A ju w rr F3 g

n=input('请输入需要计算的层数：');

rrg=[];

for y=1:n

syms (['o',num2str(y)]);

syms (['p',num2str(y)]);

rrg=[rrg,rr/(mu*m*sin(sym(['o',num2str(y)])))];

end

F3=-A*mu*m*m*(sin(sym(['o',num2str(1)]))*sin(sym(['o',num2str(1+1)]))*cos(sym(['p',num2str(1)]))*cos(sym(['p',num2str(1+1)]))+sin(sym(['o',num2str(1)]))*sin(sym(['o',num2str(1+1)]))*sin(sym(['p',num2str(1)]))*sin(sym(['p',num2str(1+1)]))+cos(sym(['o',num2str(1)]))*cos(sym(['o',num2str(1+1)])));

for  j=2

n-1)

F3=F3-A*mu*m*m*(sin(sym(['o',num2str(j)]))*sin(sym(['o',num2str(j+1)]))*cos(sym(['p',num2str(j)]))*cos(sym(['p',num2str(j+1)]))+sin(sym(['o',num2str(j)]))*sin(sym(['o',num2str(j+1)]))*sin(sym(['p',num2str(j)]))*sin(sym(['p',num2str(j+1)]))+cos(sym(['o',num2str(j)]))*cos(sym(['o',num2str(j+1)])));

end

F1=-k*sin(sym(['o',num2str(1)]))*sin(sym(['o',num2str(1)]))*cos(sym(['p',num2str(1)]))*cos(sym(['p',num2str(1)]));

F2=0.5*mu*m*m*cos(sym(['o',num2str(1)]))*cos(sym(['o',num2str(1)]));

for l=2:n

F1=F1-k*sin(sym(['o',num2str(l)]))*sin(sym(['o',num2str(l)]))*cos(sym(['p',num2str(l)]))*cos(sym(['p',num2str(l)]));

F2=F2+0.5*mu*m*m*cos(sym(['o',num2str(l)]))*cos(sym(['o',num2str(l)]));

end

F=F1+F2+F3;

for  x=1:2*n

if x==1

ju(x,1)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2))]),1)+i*w;

ju(x,2)=rrg(ceil(x/2))*diff(diff(F,sym(['p',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2))]),1);

ju(x,3)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2)+1)]),1),sym(['p',num2str(ceil(x/2))]),1);

ju(x,4)=rrg(ceil(x/2))*diff(diff(F,sym(['p',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2)+1)]),1);

continue

end

if x==2

ju(x,1)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['o',num2str(ceil(x/2))]),1);

ju(x,2)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2))]),1)-i*w;

ju(x,3)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['o',num2str(ceil(x/2)+1)]),1);

ju(x,4)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2)+1)]),1);

continue

end

if x==2*n-1

ju(x,(ceil(x/2)-2)*2+1)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2)-1)]),1),sym(['p',num2str(ceil(x/2))]),1);

ju(x,(ceil(x/2)-2)*2+2)=rrg(ceil(x/2))*diff(diff(F,sym(['p',num2str(ceil(x/2)-1)]),1),sym(['p',num2str(ceil(x/2))]),1);

ju(x,(ceil(x/2)-2)*2+3)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2))]),1)+i*w;

ju(x,(ceil(x/2)-2)*2+4)=rrg(ceil(x/2))*diff(diff(F,sym(['p',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2))]),1);

continue

end

if x==2*n

ju(x,(ceil(x/2)-2)*2+1)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2)-1)]),1),sym(['o',num2str(ceil(x/2))]),1);

ju(x,(ceil(x/2)-2)*2+2)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2)-1)]),1);

ju(x,(ceil(x/2)-2)*2+3)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['o',num2str(ceil(x/2))]),1);

ju(x,(ceil(x/2)-2)*2+4)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2))]),1)-i*w;

continue

end

xx=mod(x,2);

if(xx==1)

ju(x,(ceil(x/2)-2)*2+1)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2)-1)]),1),sym(['p',num2str(ceil(x/2))]),1);

ju(x,(ceil(x/2)-2)*2+2)=rrg(ceil(x/2))*diff(diff(F,sym(['p',num2str(ceil(x/2)-1)]),1),sym(['p',num2str(ceil(x/2))]),1);

ju(x,(ceil(x/2)-2)*2+3)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2))]),1)+i*w;

ju(x,(ceil(x/2)-2)*2+4)=rrg(ceil(x/2))*diff(diff(F,sym(['p',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2))]),1);

ju(x,(ceil(x/2)-2)*2+5)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2)+1)]),1),sym(['p',num2str(ceil(x/2))]),1);

ju(x,(ceil(x/2)-2)*2+6)=rrg(ceil(x/2))*diff(diff(F,sym(['p',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2)+1)]),1);

else

ju(x,(ceil(x/2)-2)*2+1)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2)-1)]),1),sym(['o',num2str(ceil(x/2))]),1);

ju(x,(ceil(x/2)-2)*2+2)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2)-1)]),1);

ju(x,(ceil(x/2)-2)*2+3)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['o',num2str(ceil(x/2))]),1);

ju(x,(ceil(x/2)-2)*2+4)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2))]),1)-i*w;

ju(x,(ceil(x/2)-2)*2+5)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['o',num2str(ceil(x/2)+1)]),1);

ju(x,(ceil(x/2)-2)*2+6)=rrg(ceil(x/2))*diff(diff(F,sym(['o',num2str(ceil(x/2))]),1),sym(['p',num2str(ceil(x/2)+1)]),1);

end

end

for y=1:n

syms(['o',num2str(y)]);

syms(['p',num2str(y)]);

g(2*y-1)=['o',num2str(y)];

g(2*y)=['p',num2str(y)];

d(2*y-1)=pi/2;

d(2*y)=0;

end

ji2=subs(ju,{m,k,mu,A,rr},{1.3e6,4e5,4*pi*1e-7,0.5,2.211e5});

ji3=subs(ji2,g,d);

xxx1=[0:1e8:8e11];

xxx2=zeros(8001,1);

xxx1=xxx1';

for i=1:8001

xxx=subs(ji3,w,xxx1(i));

xxx2(i)=det(xxx);

end

plot(xxx1,xxx2)