Hi, I want to plot the eigenvalue of the Orr-Sommerfeld operator at... (2024)

% p40.m - eigenvalues of Orr-Sommerfeld operator (compare p38.m)

R = 5772; clf, [ay,ax] = meshgrid([.56 .04],[.1 .5]);

for N = 40:20:100

% 2nd- and 4th-order differentiation matrices:

[D,x] = cheb(N); D2 = D^2; D2 = D2(2:N,2:N);

S = diag([0; 1 ./(1-x(2:N).^2); 0]);

D4 = (diag(1-x.^2)*D^4 - 8*diag(x)*D^3 - 12*D^2)*S;

D4 = D4(2:N,2:N);

% Orr-Sommerfeld operators A,B and generalized eigenvalues:

I = eye(N-1);

A = (D4-2*D2+I)/R - 2i*I - 1i*diag(1-x(2:N).^2)*(D2-I);

B = D2-I;

ee = eig(A,B);

i = N/20-1; subplot('position',[ax(i) ay(i) .38 .38])

plot(ee,'.','markersize',12)

grid on, axis([-.8 .2 -1 0]), axis square

title(['N = ' int2str(N) ' \lambda_{max} = ' ...

num2str(max(real(ee)),'%16.12f')]), drawnow

end

Unrecognized function or variable 'cheb'.

with the help of eigenshuffle function.

function [Vseq,Dseq] = eigenshuffle(Asequence,Bsequence)

% eigenshuffle: Consistent sorting for an eigenvalue/vector sequence

% [Vseq,Dseq] = eigenshuffle(Asequence)

% [Vseq,Dseq] = eigenshuffle(Asequence,Bsequence)

% Includes munkres.m (by gracious permission from Yi Cao)

% to choose the appropriate permutation. This greatly

% enhances the speed of eigenshuffle over my previous

% release.

% http://www.mathworks.com/matlabcentral/fileexchange/20652

% Arguments: (input)

% Asequence - an array of eigenvalue problems. If

% Asequence is a 3-d numeric array, then each

% plane of Asequence must contain a square

% matrix that will be used to call eig.

% Eig will be called on each of these matrices

% to produce a series of eigenvalues/vectors,

% one such set for each eigenvalue problem.

% Bsequence - (OPTIONAL) This allows the user to solve

% generalized eigenvalue problems of the form

% A*X = lambda*B*x

% Bsequence must be the same shape and size as

% Asequence.

% ONLY supply Bsequence when generalized eigenvalues

% will be computed.

% Arguments: (Output)

% Vseq - a 3-d array (pxpxn) of eigenvectors. Each

% plane of the array will be sorted into a

% consistent order with the other eigenvalue

% problems. The ordering chosen will be one

% that maximizes the energy of the consecutive

% eigensystems relative to each other.

% Dseq - pxn array of eigen values, sorted in order

% to be consistent with each other and with the

% eigenvectors in Vseq.

% Example:

% Efun = @(t) [1 2*t+1 t^2 t^3;2*t+1 2-t t^2 1-t^3; ...

% t^2 t^2 3-2*t t^2;t^3 1-t^3 t^2 4-3*t];

% Aseq = zeros(4,4,21);

% for i = 1:21

% Aseq(:,:,i) = Efun((i-11)/10);

% end

% [Vseq,Dseq] = eigenshuffle(Aseq);

% To see that eigenshuffle has done its work correctly,

% look at the eigenvalues in sequence, after the shuffle.

% t = (-1:.1:1)';

% [t,Dseq']

% ans =

% -1 8.4535 5 2.3447 0.20181

% -0.9 7.8121 4.7687 2.3728 0.44644

% -0.8 7.2481 4.56 2.3413 0.65054

% -0.7 6.7524 4.3648 2.2709 0.8118

% -0.6 6.3156 4.1751 2.1857 0.92364

% -0.5 5.9283 3.9855 2.1118 0.97445

% -0.4 5.5816 3.7931 2.0727 0.95254

% -0.3 5.2676 3.5976 2.0768 0.858

% -0.2 4.9791 3.3995 2.1156 0.70581

% -0.1 4.7109 3.2 2.1742 0.51494

% 0 4.4605 3 2.2391 0.30037

% 0.1 4.2302 2.8 2.2971 0.072689

% 0.2 4.0303 2.5997 2.3303 -0.16034

% 0.3 3.8817 2.4047 2.3064 -0.39272

% 0.4 3.8108 2.1464 2.2628 -0.62001

% 0.5 3.8302 1.8986 2.1111 -0.83992

% 0.6 3.9301 1.5937 1.9298 -1.0537

% 0.7 4.0927 1.2308 1.745 -1.2685

% 0.8 4.3042 0.82515 1.5729 -1.5023

% 0.9 4.5572 0.40389 1.4272 -1.7883

% 1 4.8482 -8.0012e-16 1.3273 -2.1755

if n < 2

% we can quit now, having sorted the eigenvalues

% as best as we could.

return

end

% now, treat each eigenproblem in sequence (after

% the first one.)

for i = 2:n

% compute distance between systems

V1 = Vseq(:,:,i-1);

V2 = Vseq(:,:,i);

D1 = Dseq(:,i-1);

D2 = Dseq(:,i);

dist = (1-abs(V1'*V2)).*sqrt( ...

distancematrix(real(D1),real(D2)).^2+ ...

distancematrix(imag(D1),imag(D2)).^2);

% Is there a best permutation? use munkres.

% much faster than my own mintrace, munkres

% is used by gracious permission from Yi Cao.

reorder = munkres(dist);

Vseq(:,:,i) = Vseq(:,reorder,i);

Dseq(:,i) = Dseq(reorder,i);

% also ensure the signs of each eigenvector pair

% were consistent if possible

S = squeeze(real(sum(Vseq(:,:,i-1).*Vseq(:,:,i),1))) < 0;

Vseq(:,S,i) = -Vseq(:,S,i);

end

% =================

% end mainline

% =================

% begin subfunctions

% =================

function d = distancematrix(vec1,vec2)

% simple interpoint distance matrix

[vec1,vec2] = ndgrid(vec1,vec2);

d = abs(vec1 - vec2);

function [assignment,cost] = munkres(costMat)

% MUNKRES Munkres (Hungarian) Algorithm for Linear Assignment Problem.

% [ASSIGN,COST] = munkres(COSTMAT) returns the optimal column indices,

% ASSIGN assigned to each row and the minimum COST based on the assignment

% problem represented by the COSTMAT, where the (i,j)th element represents the cost to assign the jth

% job to the ith worker.

% This is vectorized implementation of the algorithm. It is the fastest

% among all Matlab implementations of the algorithm.

% Examples

% Example 1: a 5 x 5 example

[assignment,cost] = munkres(magic(5));

disp(assignment); % 3 2 1 5 4

disp(cost); %15

% Example 2: 400 x 400 random data

n=400;

A=rand(n);

tic

[a,b]=munkres(A);

toc % about 2 seconds

% Example 3: rectangular assignment with inf costs

A=rand(10,7);

A(A>0.7)=Inf;

[a,b]=munkres(A);

% Reference:

% "Munkres' Assignment Algorithm, Modified for Rectangular Matrices",

% http://csclab.murraystate.edu/bob.pilgrim/445/munkres.html

% version 2.0 by Yi Cao at Cranfield University on 10th July 2008

assignment = zeros(1,size(costMat,1));

cost = 0;

costMat(costMat~=costMat)=Inf;

validMat = costMat<Inf;

validCol = any(validMat,1);

validRow = any(validMat,2);

nRows = sum(validRow);

nCols = sum(validCol);

n = max(nRows,nCols);

if ~n

return

end

maxv=10*max(costMat(validMat));

dMat = zeros(n) + maxv;

dMat(1:nRows,1:nCols) = costMat(validRow,validCol);

%*************************************************

% Munkres' Assignment Algorithm starts here

%*************************************************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% STEP 1: Subtract the row minimum from each row.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

minR = min(dMat,[],2);

minC = min(bsxfun(@minus, dMat, minR));

%**************************************************************************

% STEP 2: Find a zero of dMat. If there are no starred zeros in its

% column or row start the zero. Repeat for each zero

%**************************************************************************

zP = dMat == bsxfun(@plus, minC, minR);

starZ = zeros(n,1);

while any(zP(:))

[r,c]=find(zP,1);

starZ(r)=c;

zP(r,:)=false;

zP(:,c)=false;

end

while 1

%**************************************************************************

% STEP 3: Cover each column with a starred zero. If all the columns are

% covered then the matching is maximum

%**************************************************************************

if all(starZ>0)

break

end

coverColumn = false(1,n);

coverColumn(starZ(starZ>0))=true;

coverRow = false(n,1);

primeZ = zeros(n,1);

[rIdx, cIdx] = find(dMat(~coverRow,~coverColumn)==bsxfun(@plus,minR(~coverRow),minC(~coverColumn)));

while 1

%**************************************************************************

% STEP 4: Find a noncovered zero and prime it. If there is no starred

% zero in the row containing this primed zero, Go to Step 5.

% Otherwise, cover this row and uncover the column containing

% the starred zero. Continue in this manner until there are no

% uncovered zeros left. Save the smallest uncovered value and

% Go to Step 6.

%**************************************************************************

cR = find(~coverRow);

cC = find(~coverColumn);

rIdx = cR(rIdx);

cIdx = cC(cIdx);

Step = 6;

while ~isempty(cIdx)

uZr = rIdx(1);

uZc = cIdx(1);

primeZ(uZr) = uZc;

stz = starZ(uZr);

if ~stz

Step = 5;

break;

end

coverRow(uZr) = true;

coverColumn(stz) = false;

z = rIdx==uZr;

rIdx(z) = [];

cIdx(z) = [];

cR = find(~coverRow);

z = dMat(~coverRow,stz) == minR(~coverRow) + minC(stz);

rIdx = [rIdx(:);cR(z)];

cIdx = [cIdx(:);stz(ones(sum(z),1))];

end

if Step == 6

% *************************************************************************

% STEP 6: Add the minimum uncovered value to every element of each covered

% row, and subtract it from every element of each uncovered column.

% Return to Step 4 without altering any stars, primes, or covered lines.

%**************************************************************************

[minval,rIdx,cIdx]=outerplus(dMat(~coverRow,~coverColumn),minR(~coverRow),minC(~coverColumn));

minC(~coverColumn) = minC(~coverColumn) + minval;

minR(coverRow) = minR(coverRow) - minval;

else

break

end

%**************************************************************************

% STEP 5:

% Construct a series of alternating primed and starred zeros as

% follows:

% Let Z0 represent the uncovered primed zero found in Step 4.

% Let Z1 denote the starred zero in the column of Z0 (if any).

% Let Z2 denote the primed zero in the row of Z1 (there will always

% be one). Continue until the series terminates at a primed zero

% that has no starred zero in its column. Unstar each starred

% zero of the series, star each primed zero of the series, erase

% all primes and uncover every line in the matrix. Return to Step 3.

%**************************************************************************

rowZ1 = find(starZ==uZc);

starZ(uZr)=uZc;

while rowZ1>0

starZ(rowZ1)=0;

uZc = primeZ(rowZ1);

uZr = rowZ1;

rowZ1 = find(starZ==uZc);

starZ(uZr)=uZc;

end

% Cost of assignment

rowIdx = find(validRow);

colIdx = find(validCol);

starZ = starZ(1:nRows);

vIdx = starZ <= nCols;

assignment(rowIdx(vIdx)) = colIdx(starZ(vIdx));

cost = trace(costMat(assignment>0,assignment(assignment>0)));

function [minval,rIdx,cIdx]=outerplus(M,x,y)

[nx,ny]=size(M);

minval=inf;

for r=1:nx

x1=x(r);

for c=1:ny

M(r,c)=M(r,c)-(x1+y(c));

if minval>M(r,c)

minval=M(r,c);

end

[rIdx,cIdx]=find(M==minval);

This code is from SPECTRAL METHODS IN MATLAB by Trefethen. The eigenvalue problem is solved here. But for t>0, how can I calculate eigenvalues? . Dont hasitate to ask me any questien if my problem was mis understood. Thanks in advance.

Hi, I want to plot the eigenvalue of the Orr-Sommerfeld operator at... (2024)

References