请教大神关于内波的垂直模态方程的数值求解？

Lilyliu_love

clc
clear all
close all
m=201;%分层的个数，m*dz=H
H=1000;
pi=3.1415926;
N0=5.23e-3;
w=zeros(m);%垂直速度
N=zeros(m);%浮力频率
z=zeros(m-2);
f=0.0001;%科氏参量
kx=0.01;
ky=0.01;
% kx=0.01;
%ky=0.01;
dz=H/(m-1);
l=dz^2*(kx^2+ky^2);
A=zeros(m-2,m-2);
B=zeros(m-2,m-2);
C=zeros(m-2,m-2);
for i=1:m-2
   z(i)=-i*dz;
end

%N赋值
% z ≤-200m式中N(-200)一般可取为5.23×10-3秒,
%G.加勒特和W.H.蒙克在1972年发表的研究结果，在大洋中，除了高纬度海区、赤道海区和西部边界流海区外，
%上混合层以下的N(z)分布可表示为
%N(200)=5.23*10^(-3);
%考虑跃层-300m左右,取为1×10-2秒
for i=1:m-2
N(i+1)=0.01*gaussmf(-z(i)+1,[100,300])+0.005;
end
%for i=1:m
  %  N(i)=5.23e-3;
%end

% A*v=lambda*B*v
for i=1:m-2
    if i~=1&&i~=m-2
    A(i,i)=2+l*N(i+1)^2/(f*f);
    A(i,i+1)=-1;
    A(i,i-1)=-1;
    end
     if i==1
    A(i,i)=2+l*N(i+1)^2/(f*f);
    A(i,i+1)=-1;
     end
     if i==m-2
    A(i,i)=2+l*N(i+1)^2/(f*f);
    A(i,i-1)=-1;
    end
   
end
%B赋值
for i=1:m-2
    %i+2 i+1 i
     if i~=1&&i~=m-2
    B(i,i)=2+l;
    B(i,i+1)=-1;
    B(i,i-1)=-1;
    end
     if i==1
    B(i,i)=2+l;
    B(i,i+1)=-1;
     end
     if i==m-2
    B(i,i)=2+l;
    B(i,i-1)=-1;
    end
end
% A=B.*inv(C);
%
[Wv,D] = eig(A,B,'qz');
%
% % [Wv,D] = eigs(A,B,m-2);
% [Wv,D] = eig(full(A),full(B),'qz');
subplot(2,2,1 )
plot(Wv(:,1)',z)%第一模态
set(gca,'xlim',[-1,1]);
title('第一模态')
t=sqrt(D(1,1)*f*f)
% xlabel(num2str(D(1,1)))
xlabel(num2str(t))
ylabel('deepth')
subplot(2,2,2)
plot(Wv(:,2)',z)%第二模态
set(gca,'xlim',[-1,1]);
title('第二模态')
t=sqrt(D(2,2)*f*f)
xlabel(num2str(t))
% xlabel(num2str(D(2,2)))
ylabel('deepth')
subplot(2,2,3)
plot(Wv(:,3)',z)%第二模态
title('第三模态')
t=sqrt(D(3,3)*f*f)
xlabel(num2str(t))
% xlabel(num2str(D(3,3)))
ylabel('deepth')
subplot(2,2,4)
plot(N(2:m-1),z)
xlabel('N(z)')
set(gca,'xlim',[0,0.02]);
hold off;

以上是matlab程序，eig函数求解广义特征值和特征向量，作图时发现每隔一层模态就反向一次，如图：（图1 m=201也就是200层，图2为201层）

风往北吹

不太懂，默默帮顶~

wuxiangwang

不懂，默默帮顶~

1028478505

最近也在做这个遇到了一点类似问题，请问LZ问题解决了嘛