Last active
December 14, 2015 09:50
-
-
Save mathsam/5067977 to your computer and use it in GitHub Desktop.
a three layer model; calculate the linear growth rate
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
syms c epsilon k N gamma | |
A = [-(-c+0.5*(3+epsilon))*(1+k^2)+k*(1+gamma), -c+0.5*(3+epsilon), 0; ... | |
N*(-c+0.5*(1+epsilon)), -(-c+0.5*(1+epsilon))*(k^2+2*N)+k*(gamma-N*(0.5-0.5*epsilon)), N*(-c+0.5*(1+epsilon)); ... | |
0, -c*N/epsilon, c*(k^2+N/epsilon)+k*(gamma-0.5*N*(1+epsilon)/epsilon)]; | |
det_A = vpa(det(A)); | |
%% | |
clear | |
wavenumberSet = 0.01:0.02:4; | |
criticalitySet = 0.01:0.02:2; | |
c1 = zeros(length(criticalitySet),length(wavenumberSet)); | |
c2 = zeros(length(criticalitySet),length(wavenumberSet)); | |
c3 = zeros(length(criticalitySet),length(wavenumberSet)); | |
N = 0.5; | |
epsilon = 0.1; | |
for i_wavenumber = 1:length(wavenumberSet) | |
for i_criticality = 1:length(criticalitySet) | |
k = wavenumberSet(i_wavenumber); | |
gamma = criticalitySet(i_criticality); | |
poly_A(1) = (0.125*(8.*epsilon*k^4 + 8.*epsilon*k^6 + 8.*k^2*N + 8.*epsilon*k^2*N + 8.*k^4*N + 16.*epsilon*k^4*N + 8.*k^2*N^2))/epsilon; | |
poly_A(2) = (0.125*(16.*epsilon*gamma*k^3 - 16.*epsilon*k^4 - 8.*epsilon^2*k^4 + 8.*epsilon*k^5 + 24.*epsilon*gamma*k^5 - 16.*epsilon*k^6 - 8.*epsilon^2*k^6 + 8.*gamma*k*N + 8.*epsilon*gamma*k*N - 16.*k^2*N - 24.*epsilon*k^2*N - 8.*epsilon^2*k^2*N + 4.*k^3*N + 8.*epsilon*k^3*N + 4.*epsilon^2*k^3*N + 16.*gamma*k^3*N + 32.*epsilon*gamma*k^3*N - 16.*k^4*N - 40.*epsilon*k^4*N - 16.*epsilon^2*k^4*N - 4.*k^5*N - 8.*epsilon*k^5*N + 4.*epsilon^2*k^5*N + 8.*gamma*k*N^2 - 16.*k^2*N^2 - 8.*epsilon*k^2*N^2 - 12.*k^3*N^2 - 4.*epsilon*k^3*N^2))/epsilon; | |
poly_A(3) = (0.125*(8.*epsilon*gamma^2*k^2 - 28.*epsilon*gamma*k^3 - 12.*epsilon^2*gamma*k^3 + 6.*epsilon*k^4 + 8.*epsilon^2*k^4 + 2.*epsilon^3*k^4 + 16.*epsilon*gamma*k^4 + 24.*epsilon*gamma^2*k^4 - 4.*epsilon*k^5 - 4.*epsilon^2*k^5 - 32.*epsilon*gamma*k^5 - 16.*epsilon^2*gamma*k^5 + 6.*epsilon*k^6 + 8.*epsilon^2*k^6 + 2.*epsilon^3*k^6 - 12.*gamma*k*N - 20.*epsilon*gamma*k*N - 8.*epsilon^2*gamma*k*N + 6.*k^2*N + 14.*epsilon*k^2*N + 10.*epsilon^2*k^2*N + 2.*epsilon^3*k^2*N + 4.*gamma*k^2*N + 8.*epsilon*gamma*k^2*N + 4.*epsilon^2*gamma*k^2*N + 8.*gamma^2*k^2*N + 16.*epsilon*gamma^2*k^2*N + 4.*k^3*N + 6.*epsilon*k^3*N - 8.*epsilon^2*k^3*N - 2.*epsilon^3*k^3*N - 16.*gamma*k^3*N - 48.*epsilon*gamma*k^3*N - 24.*epsilon^2*gamma*k^3*N + 2.*k^4*N + 12.*epsilon*k^4*N + 22.*epsilon^2*k^4*N + 4.*epsilon^3*k^4*N - 8.*gamma*k^4*N - 16.*epsilon*gamma*k^4*N + 8.*epsilon^2*gamma*k^4*N + 8.*k^5*N + 18.*epsilon*k^5*N - 2.*epsilon^3*k^5*N + 10.*k*N^2 + 4.*epsilon*k*N^2 + 2.*epsilon^2*k*N^2 - 4.*gamma*k*N^2 - 4.*epsilon*gamma*k*N^2 - 4.*k^2*N^2 + 4.*epsilon*k^2*N^2 - 12.*gamma*k^2*N^2 - 4.*epsilon*gamma*k^2*N^2 + 22.*k^3*N^2 + 20.*epsilon*k^3*N^2 + 6.*epsilon^2*k^3*N^2 + 2.*k^4*N^2 - 2.*epsilon^2*k^4*N^2))/epsilon; | |
poly_A(4) = (0.125*(-12.*epsilon*gamma^2*k^2 - 4.*epsilon^2*gamma^2*k^2 + 6.*epsilon*gamma*k^3 + 8.*epsilon^2*gamma*k^3 + 2.*epsilon^3*gamma*k^3 + 8.*epsilon*gamma^2*k^3 + 8.*epsilon*gamma^3*k^3 - 4.*epsilon*gamma*k^4 - 4.*epsilon^2*gamma*k^4 - 16.*epsilon*gamma^2*k^4 - 8.*epsilon^2*gamma^2*k^4 + 6.*epsilon*gamma*k^5 + 8.*epsilon^2*gamma*k^5 + 2.*epsilon^3*gamma*k^5 + 6.*epsilon*gamma*k*N + 8.*epsilon^2*gamma*k*N + 2.*epsilon^3*gamma*k*N + 6.*gamma*k^2*N + 6.*epsilon*gamma*k^2*N - 10.*epsilon^2*gamma*k^2*N - 2.*epsilon^3*gamma*k^2*N - 8.*epsilon*gamma^2*k^2*N - 8.*epsilon^2*gamma^2*k^2*N - 3.*k^3*N - 7.*epsilon*k^3*N - 5.*epsilon^2*k^3*N - 1.*epsilon^3*k^3*N - 4.*gamma*k^3*N + 4.*epsilon*gamma*k^3*N + 20.*epsilon^2*gamma*k^3*N + 4.*epsilon^3*gamma*k^3*N - 4.*gamma^2*k^3*N - 8.*epsilon*gamma^2*k^3*N + 4.*epsilon^2*gamma^2*k^3*N + 2.*k^4*N + 4.*epsilon*k^4*N + 2.*epsilon^2*k^4*N + 8.*gamma*k^4*N + 18.*epsilon*gamma*k^4*N - 2.*epsilon^3*gamma*k^4*N - 3.*k^5*N - 7.*epsilon*k^5*N - 5.*epsilon^2*k^5*N - 1.*epsilon^3*k^5*N - 3.*k*N^2 - 7.*epsilon*k*N^2 - 5.*epsilon^2*k*N^2 - 1.*epsilon^3*k*N^2 + k^2*N^2 + 7.*epsilon*k^2*N^2 + 7.*epsilon^2*k^2*N^2 + epsilon^3*k^2*N^2 + 4.*gamma*k^2*N^2 + 8.*epsilon*gamma*k^2*N^2 + 4.*epsilon^2*gamma*k^2*N^2 - 4.*k^3*N^2 - 14.*epsilon*k^3*N^2 - 12.*epsilon^2*k^3*N^2 - 2.*epsilon^3*k^3*N^2 + 2.*gamma*k^3*N^2 - 2.*epsilon^2*gamma*k^3*N^2 - 3.*k^4*N^2 - 1.*epsilon*k^4*N^2 + 3.*epsilon^2*k^4*N^2 + epsilon^3*k^4*N^2))/epsilon; | |
waveSpeed = roots(poly_A); | |
c1(i_criticality,i_wavenumber) = waveSpeed(1); | |
c2(i_criticality,i_wavenumber) = waveSpeed(2); | |
c3(i_criticality,i_wavenumber) = waveSpeed(3); | |
end | |
end | |
%% plot | |
wavenumber2d = repmat(wavenumberSet(:)',length(criticalitySet),1); | |
tmp1 = imag(c1.*wavenumber2d); | |
tmp2 = imag(c2.*wavenumber2d); | |
tmp3 = imag(c3.*wavenumber2d); | |
tmpindex = tmp1(:)>0; | |
growthRate1 = zeros(size(tmp1)); | |
growthRate2 = zeros(size(tmp1)); | |
growthRate1(tmpindex) = tmp1(tmpindex); | |
tmpindex = tmp2(:)>0; | |
growthRate1(tmpindex) = tmp2(tmpindex); | |
tmpindex = tmp2(:)<0; | |
growthRate2(tmpindex) = tmp2(tmpindex); | |
tmpindex = tmp3(:)<0; | |
growthRate2(tmpindex) = tmp3(tmpindex); | |
figure | |
contourf(wavenumberSet,criticalitySet,growthRate1,30) | |
colorbar | |
set(gca,'YDir','reverse') | |
xlabel('wavenumber','fontsize',16) | |
ylabel('1/Sc','fontsize',16) | |
set(gca,'fontsize',16) | |
xlim([0 2]) | |
figure | |
contourf(wavenumberSet,criticalitySet,growthRate2,30) | |
colorbar | |
set(gca,'YDir','reverse') | |
xlabel('wavenumber','fontsize',16) | |
ylabel('1/Sc','fontsize',16) | |
set(gca,'fontsize',16) | |
xlim([0 2]) | |
%ylim([0.2 0.5]) | |
%% calculate the eigen vectors correspond to phase speed (growthRate>0) | |
wavenumber2d = repmat(wavenumberSet(:)',length(criticalitySet),1); | |
tmp1 = c1.*wavenumber2d; | |
tmp2 = c2.*wavenumber2d; | |
tmpindex = imag(tmp1(:))>0; | |
mod1 = zeros(size(tmp1)); | |
mod1(tmpindex) = tmp1(tmpindex); | |
tmpindex = imag(tmp2(:))>0; | |
mod1(tmpindex) = tmp2(tmpindex); | |
wavenumberSet = 0.01:0.02:4; | |
criticalitySet = 0.01:0.02:2; | |
N = 0.5; | |
epsilon = 0.1; | |
phi = zeros(size(c1,1),size(c1,2),9); | |
for i_wavenumber = 1:length(wavenumberSet) | |
for i_criticality = 1:length(criticalitySet) | |
c = mod1(i_criticality,i_wavenumber); | |
k = wavenumberSet(i_wavenumber); | |
gamma = criticalitySet(i_criticality); | |
A = [-(-c+0.5*(3+epsilon))*(1+k^2)+k*(1+gamma), -c+0.5*(3+epsilon), 0; ... | |
N*(-c+0.5*(1+epsilon)), -(-c+0.5*(1+epsilon))*(k^2+2*N)+k*(gamma-N*(0.5-0.5*epsilon)), N*(-c+0.5*(1+epsilon)); ... | |
0, -c*N/epsilon, c*(k^2+N/epsilon)+k*(gamma-0.5*N*(1+epsilon)/epsilon)]; | |
[V D] = eig(A); | |
phi(i_criticality,i_wavenumber,1:3) = V(1:3)/V(3); | |
phi(i_criticality,i_wavenumber,4:6) = V(4:6)/V(6); | |
phi(i_criticality,i_wavenumber,7:9) = V(7:9)/V(9); | |
end | |
end | |
%% | |
subplot(3,3,1) | |
imagesc(real(phi(:,:,1))) | |
subplot(3,3,2) | |
imagesc(real(phi(:,:,2))) | |
subplot(3,3,3) | |
imagesc(real(phi(:,:,3))) | |
subplot(3,3,4) | |
imagesc(real(phi(:,:,4))) | |
subplot(3,3,5) | |
imagesc(real(phi(:,:,5))) | |
subplot(3,3,6) | |
imagesc(real(phi(:,:,6))) | |
subplot(3,3,7) | |
imagesc(real(phi(:,:,7))) | |
subplot(3,3,8) | |
imagesc(real(phi(:,:,8))) | |
subplot(3,3,9) | |
imagesc(real(phi(:,:,9))) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment