First commit of calculating the eigenvectors and values manually.

Author: AndreyAPopov

src/+otp/+qg/QuasiGeostrophicProblem.m

@@ -117,7 +117,7 @@ function onSettingsChanged(obj)
             n = [nx, ny];
 
             hx = 1/(nx + 1);
-            hy = 1/(ny + 1);
+            hy = 2/(ny + 1);
 
             xdomain = [0, 1];
             ydomain = [0, 2];
@@ -152,23 +152,22 @@ function onSettingsChanged(obj)
             % for a detailed method, the "Eigenvalue Method" which makes this
             % particularly efficient
 
-            nfLx = -full(Lx);
-            nfLy = -full(Ly);
-            [P1, Lambda] = eig(nfLx);
-            [P2, D] = eig(nfLy);
+            %nfLx = -full(Lx);
+            %nfLy = -full(Ly);
+            %[P1, Lambda] = eig(nfLx);
+            %[P2, D] = eig(nfLy);
 
             % We can represent the eigenvalues as
-            %dLambda = (4/(hx^2) * (sin(pi*(1:nx)/(2*(nx + 1))).^2)).';
-            %dD = (4/(hy^2) * (sin(pi*(1:ny)/(2*(ny + 1))).^2)).';
-            %L12 = 1./(dLambda + dD.');
-            %P1 = sqrt(2/(nx + 1))*sin((1:nx).'*(1:nx)*pi/(nx + 1));
-            %P2 = sqrt(2/(ny + 1))*sin((1:ny).'*(1:ny)*pi/(ny + 1));
+            dLambda = (4/(hx^2) * (sin(pi*(1:nx)/(2*(nx + 1))).^2)).';
+            dD = (4/(hy^2) * (sin(pi*(1:ny)/(2*(ny + 1))).^2)).';
+            L12 = 1./(dLambda + dD.');
+            P1 = sqrt(2/(nx + 1))*sin((1:nx).'*(1:nx)*pi/(nx + 1));
+            P2 = sqrt(2/(ny + 1))*sin((1:ny).'*(1:ny)*pi/(ny + 1));
 
-            
-            
-            L12 = 1./(diag(Lambda) + diag(D).');
-            P1T = P1.';
-            P2T = P2.';
+
+            %L12 = 1./(diag(Lambda) + diag(D).');
+            %P1T = P1.';
+            %P2T = P2.';
 
             ys = linspace(ydomain(1), ydomain(end), ny + 2);
             ys = ys(2:end-1);
@@ -179,7 +178,7 @@ function onSettingsChanged(obj)
             F = sin(pi*(ymat.' - 1));
 
             obj.Rhs = otp.Rhs(@(t, psi) ...
-                otp.qg.f(psi, Lx, Ly, P1, P1T, P2, P2T, L12, Dx, DyT, F, Re, Ro), ...
+                otp.qg.f(psi, Lx, Ly, P1, P2, L12, Dx, DyT, F, Re, Ro), ...
                 ...
                 otp.Rhs.FieldNames.JacobianVectorProduct, @(t, psi, u) ...
                 otp.qg.jvp(psi, u, L, RdnL, PdnL, Ddx, Ddy, Re, Ro), ...

src/+otp/+qg/f.m

@@ -1,4 +1,4 @@
-function dpsit = f(psi, Lx, Ly, P1, P1T, P2, P2T, L12, Dx, DyT, F, Re, Ro)
+function dpsit = f(psi, Lx, Ly, P1, P2, L12, Dx, DyT, F, Re, Ro)
 
 [nx, ny] = size(L12);
 
@@ -24,7 +24,7 @@
 dqtmq = mJ + (1/Ro)*(dpsix) + (1/Ro)*F;
 
 % solve the sylvester equation
-nLidqtmq = P1*(L12.*(P1T*dqtmq*P2))*P2T;
+nLidqtmq = P1*(L12.*(P1*dqtmq*P2))*P2;
 
 % solve into stream form of the rhs
 dpsit = reshape(nLidqtmq - (1/Re)*(q), nx*ny, 1);