# 1

Apr 20th, 2021
647
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. clc; close all
2.
3. % [U, S, V] = svd(Asnap, 'econ');
4. % plot(diag(S)/sum(diag(S)), 'ro')
5. % % based on this plot, we set r = 10
6.
7.
8. %% DMD
9. X1 = Asnap(:, 1:end-1);
10. X2 = Asnap(:, 2:end);
11. r = 10;
12.
13. [U, S, V] = svd(X1, 'econ');
14.
15. % rank-r truncation
16. Ur = U(:, 1:r);
17. Sr = S(1:r, 1:r);
18. Vr = V(:, 1:r);
19.
20. Atilde = Ur'*X2*Vr/Sr;
21.
22. % eigenvalues and eigenvectors
23. [W, D] = eig(Atilde);
24.
25. % DMD modes
26. Phi = X2*Vr/Sr*W;
27.
28. % eigenvalues
29. lamdba = diag(D);
30.
31. % frequencies
32. omega = log(lamdba)/(iconv*Deltat);
33.
34. % initial condition
35. x1 = Asnap(:, 1);
36. b = Phi\x1;
37.
38. %% Reconstruction
39. time_vector = (0:iconv*Deltat:nt*Deltat);
40. time_dynamics = zeros(r, size(Asnap, 2));
41. for iter = 1:size(Asnap, 2)
42.     time_dynamics(:, iter) = (b.*exp(omega*time_vector(iter)));
43. end
44. Asnap_dmd = Phi*time_dynamics;
45.
46. % get vorticity vector at last time step
47. vorticity_dmd = real(Asnap_dmd(:, 100));
48.
49. % reshape for plotting
50. vorticity_dmd = reshape(vorticity_dmd, ncx-1, ncy-1);
51. vorticity_dmd = vorticity_dmd';
52.
53. % plot
54. fz=figure;
55. az=axes;
56. maxval=max(max(vorticity_dmd));
57. minval=min(min(vorticity_dmd));
58. vals=linspace(minval,maxval,100);
59. contour(X,Y,vorticity_dmd,vals);
60. set(az,'DataAspectRatio',[1 1 1],'XLim',[0 lx],'YLim',[0 ly]);
61. title(sprintf(' Vorticity (DMD) (%g...%g); Re_b=%g',...
62.                 minval,maxval,Re));
63. colorbar;
64.
RAW Paste Data