请各位高手指点我的理解那个地方有错?
本帖最后由 free518 于 2012-12-7 06:49 编辑我在看一篇文献时想实现文章中的图,结果matlab程序运行的结果怎么也不出现人家的结果,不知道哪个地方有问题?请各位高手指点路?谢谢!可能是我理解这个文献有错误。我的理解是：（1）当参数s=0，α=0时，有奇数阶矩阵，矩阵中的，
An中的n=1,2,3......
约定qn是根的正实部分. 当α=n=0时，上角标的*代表复共轭。想要求这个矩阵的特征值，再令a[i]为第i个最大特征值的倒数。i=1,2,3,4.（2）当参数s=0，α=1/2时，有偶数阶矩阵，想要求这个矩阵的特征值，再令b[i]为第i个最大特征值的倒数。i=1,2,3,4.后想把上述(1)和(2)画在一个图上，画x轴为k, y轴为a/g和b/g的图，见文献
（Linear theory of Faraday instability in viscous liquids的图一）其他的参数为：ν=1.02×10-4m2s-1,σ=67.6×10-3Nm-1, h=2×10-3m, ω/(2π)=60HZ , ρ=1000kg/m3
我写的MatLab程序:
%% ===============================
clear all
clc
format short
Jn=31;%矩阵的阶数
omega=2*pi*60;
s=0;
alpha=0.0;
beta=0.5;
g=9.8*10^(3);
xima=67.6*10^(-3)*10^(-3);
h=2*10^(-3)*10^(3);
mu=1.02*10^(-4)*10^(6);
ro=1*10^3*10^(-9);
nn=floor(Jn/2);
i=sqrt(-1);
X_Start=0;X_end=3;Y_Start=0;Y_end=150;
for k=0.01:0.01:X_end
for n=0:nn
if n==0
AN(n+1)=(2/k)*(g*k+xima/ro*k^3);%AN(1)=A_0
else
q(n+1)=(sqrt(k^2+(s+i*(alpha+n)*omega)/mu));%q(1)=q_0,q(2)=q_1,
C(n+1)=q(n+1)*(q(n+1)^4+2*q(n+1)^2*k^2+5*k^4);
D(n+1)=k*(q(n+1)^4+6*q(n+1)^2*k^2+k^4);
AN(n+1)=(2/k)*(g*k+k^3*xima/ro-mu^2*((4*q(n+1)*k^2*(q(n+1)^2+k^2)...
-C(n+1)*cosh(q(n+1)*h)*cosh(k*h)+D(n+1)*sinh(q(n+1)*h)*sinh(k*h))...
/(q(n+1)*cosh(q(n+1)*h)*sinh(k*h)-k*sinh(q(n+1)*h)*cosh(k*h)))); %AN(2)=A_1,AN(3)=A_2
end
end
A=zeros(Jn);
for ij=1:Jn
if ij==1
A(ij,ij+1)=1/conj(AN(nn+1));
elseif ij==Jn
A(ij,ij-1)=1/AN(nn+1);
elseif (nn-ij)>=0
A(ij,ij-1)=1/conj(AN(nn-ij+2));
A(ij,ij+1)=1/conj(AN(nn-ij+2));
else
A(ij,ij-1)=1/AN(ij-nn);
A(ij,ij+1)=1/AN(ij-nn);
end
end
ab=sort(real(eig(A)),'descend');
if ab~=zeros(length(ab),1)
a=1./ab(ab~=0)/g;
plot(k,a(1),'k.',k,a(2),'k.',...
k,a(3),'k.',k,a(4),'k.',...
'LineWidth',2);
hold on
end
for n=0:nn-1
q(n+1)=(sqrt(k^2+(s+i*(beta+n)*omega)/mu));%q(1)=q_0,q(2)=q_1,
C(n+1)=q(n+1)*(q(n+1)^4+2*q(n+1)^2*k^2+5*k^4);
D(n+1)=k*(q(n+1)^4+6*q(n+1)^2*k^2+k^4);
AN(n+1)=(2/k)*(g*k+k^3*xima/ro-mu^2*((4*q(n+1)*k^2*(q(n+1)^2+k^2)...
-C(n+1)*cosh(q(n+1)*h)*cosh(k*h)+D(n+1)*sinh(q(n+1)*h)*sinh(k*h))...
/(q(n+1)*cosh(q(n+1)*h)*sinh(k*h)-k*sinh(q(n+1)*h)*cosh(k*h)))); %A(1)=A_0,A(2)=A_1
end
B=zeros(Jn-1);
for ij=1:Jn-1
if ij==1
B(ij,ij+1)=1/conj(AN(nn));
elseif ij==Jn-1
B(ij,ij-1)=1/AN(nn);
elseif (nn-ij)>=0
B(ij,ij-1)=1/conj(AN(nn-ij+1));
B(ij,ij+1)=1/conj(AN(nn-ij+1));
else
B(ij,ij-1)=1/AN(ij-nn);
B(ij,ij+1)=1/AN(ij-nn);
end
end
abab=sort(real(eig(B)),'descend');
if abab~=zeros(length(abab),1)
b=1./abab(abab~=0)/g;
plot(k,b(1),'m.',k,b(2),'m.',...
k,b(3),'m.',k,b(4),'m.',...
'LineWidth',2);
hold on
end
end
xlabel('k/{mm^{-1}}');
ylabel('a/g');
axis();
本帖最后由 free518 于 2012-12-6 23:41 编辑
文献上的图
我画出来的图
不知道哪个地方错了?请知道的高手指点,谢谢!
怎么到目前为止一个人也没有留言，难道这方面没有人考虑过吗！？
页:
[1]