请各位高手指点我的理解那个地方有错?

free518

我在看一篇文献时想实现文章中的图,结果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

Linear Theory of Faraday Instability in Viscous Liquids.pdf (309.6 KB, 下载次数: 7)

2012-12-6 23:26 上传

点击文件名下载附件
Linear theory of Faraday instability in viscous liquids

我写的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([X_Start,X_end,Y_Start,Y_end]);

free518

文献上的图

我画出来的图

不知道哪个地方错了?请知道的高手指点,谢谢!

free518

怎么到目前为止一个人也没有留言，难道这方面没有人考虑过吗！？