diff --git a/2DLinearTimeVarying/2DTimeVarying.m b/2DLinearTimeVarying/2DTimeVarying.m new file mode 100644 index 0000000..fdb01ec --- /dev/null +++ b/2DLinearTimeVarying/2DTimeVarying.m @@ -0,0 +1,196 @@ + +%% clean workspace +clear all +close all +clc + +%% define dynamics +epsilon = 1e-1; +dyn = @(t,x) ([0, 1+epsilon*t; -(1+epsilon*t),0])*x; + +% generate data +dt = 1e-1; +tspan = 0:dt:10; +x0 = [1;0]; +[tq,xq] = ode45(dyn, tspan, x0); +% extract snapshot pairs +xq = xq'; tq = tq'; +x = xq(:,1:end-1); y = xq(:,2:end); t = tq(2:end); + +% true dynamics, eigenvalues and modes +[n, m] = size(x); +A = zeros(2,2,m); +evals = zeros(2,m); +for k = 1:m + A(:,:,k) = [0, 1+epsilon*t(k); -(1+epsilon*t(k)),0]; % continuous time dynamics + evals(:,k) = eig(A(:,:,k)); +end + +%% visualize snapshots +figure +plot(tq,xq(1,:),'k-',tq,xq(2,:),'k--','MarkerSize',12,'LineWidth',1.5) +xlabel('Time $t$','Interpreter','latex') +ylabel('$x_1(t), x_2(t)$','Interpreter','latex') +fl = legend('$x_1$', '$x_2$'); +set(fl,'Interpreter','latex','Box','off','FontSize',20,'Location','best'); +box on +set(gca,'FontSize',20,'LineWidth',1) + +%% run and compare different DMD algorithms +w = 10; + +% batch DMD +AbatchDMD = zeros(2,2,m); +evalsbatchDMD = zeros(2,m); +for k = 1:w + AbatchDMD(:,:,k) = expm(A(:,:,k)*dt); +end +tic +% start at w+1 +for k = w+1:m + AbatchDMD(:,:,k) = y(:,1:k)*pinv(x(:,1:k)); +end +batchDMDelapsedTime = toc +for k = 1:m + evalA = eig(AbatchDMD(:,:,k)); + evalsbatchDMD(:,k) = log(evalA)/dt; +end + +% online DMD, rho = 1 +AonlineDMD = zeros(2,2,m); +evalsonlineDMD = zeros(2,m); +for k = 1:w + AonlineDMD(:,:,k) = expm(A(:,:,k)*dt); +end +% initial condition +odmd = OnlineDMD(2,1); +odmd.initialize(x(:,1:w),y(:,1:w)); +tic +% start at w+1 +for k = w+1:m + odmd.update(x(:,k),y(:,k)); + AonlineDMD(:,:,k) = odmd.A; +end +onlineDMDelapsedTime = toc +for k = 1:m + evalA = eig(AonlineDMD(:,:,k)); + evalsonlineDMD(:,k) = log(evalA)/dt; +end +evalsonlineDMD1 = evalsonlineDMD; + +% online DMD, rho = 0.95 +AonlineDMD = zeros(2,2,m); +evalsonlineDMD = zeros(2,m); +for k = 1:w + AonlineDMD(:,:,k) = expm(A(:,:,k)*dt); +end +% initial condition +odmd = OnlineDMD(2,0.95); +odmd.initialize(x(:,1:w),y(:,1:w)); +tic +% start at w+1 +for k = w+1:m + odmd.update(x(:,k),y(:,k)); + AonlineDMD(:,:,k) = odmd.A; +end +onlineDMDelapsedTime = toc +for k = 1:m + evalA = eig(AonlineDMD(:,:,k)); + evalsonlineDMD(:,k) = log(evalA)/dt; +end +evalsonlineDMD2 = evalsonlineDMD; + +% online DMD, rho = 0.8 +AonlineDMD = zeros(2,2,m); +evalsonlineDMD = zeros(2,m); +for k = 1:w + AonlineDMD(:,:,k) = expm(A(:,:,k)*dt); +end +% initial condition +odmd = OnlineDMD(2,0.8); +odmd.initialize(x(:,1:w),y(:,1:w)); +tic +% start at w+1 +for k = w+1:m + odmd.update(x(:,k),y(:,k)); + AonlineDMD(:,:,k) = odmd.A; +end +onlineDMDelapsedTime = toc +for k = 1:m + evalA = eig(AonlineDMD(:,:,k)); + evalsonlineDMD(:,k) = log(evalA)/dt; +end +evalsonlineDMD3 = evalsonlineDMD; + +% streaming DMD +evalsstreamingDMD = zeros(2,m); +% initial +evalsstreamingDMD(:,1) = evals(:,1); +% streaming DMD +max_rank = 2; +sdmd = StreamingDMD(max_rank); +for k = 2:m + sdmd = sdmd.update(x(:,k), y(:,k)); + [~, evalA] = sdmd.compute_modes(); + evalsstreamingDMD(:,k) = log(evalA)/dt; +end + +% mini-batch DMD +AminibatchDMD = zeros(2,2,m); +evalsminibatchDMD = zeros(2,m); +for k = 1:w + AminibatchDMD(:,:,k) = expm(A(:,:,k)*dt); +end +tic +% start at w+1 +for k = w+1:m + AminibatchDMD(:,:,k) = y(:,k-w+1:k)*pinv(x(:,k-w+1:k)); +end +minibatchDMDelapsedTime = toc +for k = 1:m + evalA = eig(AminibatchDMD(:,:,k)); + evalsminibatchDMD(:,k) = log(evalA)/dt; +end + +% window DMD +AwindowDMD = zeros(2,2,m); +evalswindowDMD = zeros(2,m); +for k = 1:w + AwindowDMD(:,:,k) = expm(A(:,:,k)*dt); +end +% initialization +wdmd = WindowDMD(2,w,1); +wdmd.initialize(x(:,1:w), y(:,1:w)); +tic +% start at w+1 +for k = w+1:m + wdmd.update(x(:,k), y(:,k)); + AwindowDMD(:,:,k) = wdmd.A; +end +windowDMDelapsedTime = toc +for k = 1:m + evalA = eig(AwindowDMD(:,:,k)); + evalswindowDMD(:,k) = log(evalA)/dt; +end + +%% Compare various DMD algorithm +index = w+1:length(t); + +% plot Imaginary part + +figure, hold on +plot(t,imag(evals(1,:)),'k-','LineWidth',1) +plot(t(1,index),imag(evalsbatchDMD(1,index)),'s-','LineWidth',1,'MarkerSize',12,'MarkerIndices',1:10:length(index)) +plot(t(1,index),imag(evalsminibatchDMD(1,index)),'o-','LineWidth',1,'MarkerSize',12,'MarkerIndices',1:10:length(index)) +plot(t(1,index),imag(evalsstreamingDMD(1,index)),'x-','LineWidth',1,'MarkerSize',12,'MarkerIndices',7:10:length(index)) +plot(t(1,index),imag(evalsonlineDMD1(1,index)),'>-','LineWidth',1,'MarkerSize',12,'MarkerIndices',4:10:length(index)) +plot(t(1,index),imag(evalsonlineDMD2(1,index)),'*-','LineWidth',1,'MarkerSize',12,'MarkerIndices',1:10:length(index)) +plot(t(1,index),imag(evalsonlineDMD3(1,index)),'d-','LineWidth',1,'MarkerSize',12,'MarkerIndices',4:10:length(index)) +plot(t(1,index),imag(evalswindowDMD(1,index)),'+-','LineWidth',1,'MarkerSize',12,'MarkerIndices',7:10:length(index)) +xlabel('Time $t$','Interpreter','latex'), ylabel('Im$(\lambda)$','Interpreter','latex') +fl = legend('true','batch','mini-batch','streaming','online, $\rho=1$','online, $\rho=0.95$','online, $\rho=0.8$','window'); +set(fl,'Interpreter','latex','Location','northwest','FontSize',20,'Box','off'); +xlim([0,10]), ylim([1,2]) + +box on +set(gca,'FontSize',20,'LineWidth',1) diff --git a/2DLinearTimeVarying/OnlineDMD.m b/2DLinearTimeVarying/OnlineDMD.m new file mode 100644 index 0000000..0959042 --- /dev/null +++ b/2DLinearTimeVarying/OnlineDMD.m @@ -0,0 +1,129 @@ +% OnlineDMD is a class that implements online dynamic mode decomposition +% The time complexity (multiply-add operation for one iteration) is O(4n^2) +% , and space complexity is O(2n^2), where n is the state dimension. +% +% Algorithm description: +% At time step k, define two matrix X(k) = [x(1),x(2),...,x(k)], +% Y(k) = [y(1),y(2),...,y(k)], that contain all the past snapshot +% pairs, where x(k), y(k) are the n dimensional state vector, +% y(k) = f(x(k)) is the image of x(k), f() is the dynamics. +% Here, if the (discrete-time) dynamics are given by z(k) = f(z(k-1)) +% , then x(k), y(k) should be measurements corresponding to +% consecutive states z(k-1) and z(k). +% At time step k+1, we need to include new snapshot pair x(k+1), y(k+1) +% We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively +% by efficient rank-1 updating online DMD algorithm. +% An exponential weighting factor can be used to place more weight on +% recent data. +% +% Usage: +% odmd = OnlineDMD(n,weighting) +% odmd.initialize(Xq,Yq) +% odmd.initilizeghost() +% odmd.update(x,y) +% [evals, modes] = odmd.computemodes() +% +% properties: +% n: state dimension +% weighting: weighting factor in (0,1] +% timestep: number of snapshot pairs processed +% A: DMD matrix, size n by n +% P: matrix that contains information about past snapshots, size n by n +% +% methods: +% initialize(Xq, Yq), initialize online DMD algorithm with q snapshot +% pairs stored in (Xq, Yq) +% initializeghost(), initialize online DMD algorithm with epsilon +% small (1e-15) ghost snapshot pairs before t=0 +% update(x,y), update when new snapshot pair (x,y) becomes available +% Here, if the (discrete-time) dynamics are given by +% z(k) = f(z(k-1)), then (x,y) should be measurements +% correponding to consecutive states z(k-1) and z(k). +% computemodes(), compute and return DMD eigenvalues and DMD modes +% +% Authors: +% Hao Zhang +% Clarence W. Rowley +% +% Reference: +% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, +% ``Online Dynamic Mode Decomposition for Time-varying Systems", +% in production, 2017. Available on arXiv. +% +% Created: +% April 2017. +% +% To look up the documentation in the command window, type help OnlineDMD + +classdef OnlineDMD < handle + properties + n = 0; % state dimension + weighting = 1; % weighting factor + timestep = 0; % number of snapshots processed + A; % DMD matrix + P; % matrix that contains information about past snapshots + end + + methods + function obj = OnlineDMD(n,weighting) + % Creat an object for online DMD + % Usage: odmd = OnlineDMD(n,weighting) + if nargin == 2 + obj.n = n; + obj.weighting = weighting; + obj.A = zeros(n,n); + obj.P = zeros(n,n); + end + end + + function initialize(obj, Xq, Yq) + % Initialize OnlineDMD with q snapshot pairs stored in (Xq, Yq) + % Usage: odmd.initialize(Xq,Yq) + q = length(Xq(1,:)); + if(obj.timestep == 0 && rank(Xq) == obj.n) + weight = (sqrt(obj.weighting)).^(q-1:-1:0); + Xq = Xq.*weight; + Yq = Yq.*weight; + obj.A = Yq*pinv(Xq); + obj.P = inv(Xq*Xq')/obj.weighting; + end + obj.timestep = obj.timestep + q; + end + + function initializeghost(obj) + % Initialize online DMD with epsilon small (1e-15) ghost + % snapshot pairs before t=0 + % Usage: odmd.initilizeghost() + epsilon = 1e-15; + obj.A = randn(obj.n, obj.n); + obj.P = (1/epsilon)*eye(obj.n); + end + + function update(obj, x, y) + % Update the DMD computation with a new pair of snapshots (x,y) + % Here, if the (discrete-time) dynamics are given by z(k) = + % f(z(k-1)), then (x,y) should be measurements correponding to + % consecutive states z(k-1) and z(k). + % Usage: odmd.update(x, y) + + % compute P*x matrix vector product beforehand + Px = obj.P*x; + % Compute gamma + gamma = 1/(1+x'*Px); + % Update A + obj.A = obj.A + (gamma*(y-obj.A*x))*Px'; + % Update P, group Px*Px' to ensure positive definite + obj.P = (obj.P - gamma*(Px*Px'))/obj.weighting; + % ensure P is SPD by taking its symmetric part + obj.P = (obj.P+(obj.P)')/2; + % time step + 1 + obj.timestep = obj.timestep + 1; + end + + function [evals, modes] = computemodes(obj) + % Compute DMD eigenvalues and DMD modes at current time + % Usage: [evals, modes] = odmd.modes() + [modes, evals] = eig(obj.A, 'vector'); + end + end +end \ No newline at end of file diff --git a/2DLinearTimeVarying/StreamingDMD.m b/2DLinearTimeVarying/StreamingDMD.m new file mode 100644 index 0000000..679f264 --- /dev/null +++ b/2DLinearTimeVarying/StreamingDMD.m @@ -0,0 +1,174 @@ +%STREAMINGDMD incrementally updated Dynamic Mode Decomposition +% Usage: +% sdmd = StreamingDMD +% sdmd = StreamingDMD(r) +% +% r: the maximum dimension of the subspace on which DMD modes are +% computed, using a POD-truncation version of the algorithm. +% If r=0 (the default), use the direct algorithm, equivalent +% to a batch-processed DMD. +% +% A StreamingDMD object has two methods: +% update adds a pair of snapshots to the dataset +% compute_modes computes the DMD modes and eigenvalues +% +% For details, see +% M.S. Hemati, M.O. Williams, C.W. Rowley, +% ``Dynamic Mode Decomposition for Large and Streaming Datasets,'' +% submitted to Physics of Fluids (2014). +% +% Tested using Matlab 8.1.0.604 (R2013a) for Mac OSX. +% +% Authors: +% Maziar S. Hemati +% Matthew O. Williams +% Clarence W. Rowley +% +% Created: +% August 2014. +% +% See also: StreamingDMD/UPDATE, StreamingDMD/COMPUTE_MODES, sdmd_run +% +% Reference page in Help browser: +% doc StreamingDMD + +classdef StreamingDMD < handle + properties + count = 0; % number of pairs of snapshots processed + max_rank = 0; % maximum allowable rank of DMD operator + Qx; % orthonormal basis for the image of X (see Eq. (3)) + Qy; % orthonormal basis for the image of Y (see Eq. (3)) + A; % the matrix $A$ needed to compute $\tilde{K}$ (see Eq. (3)) + Gx; % the matrix $G_X$ needed to compute $\tilde{K}$ (see Eq. (3)) + Gy; % the matrix $G_Y$ needed to project onto a truncated POD basis + end + + methods + function obj = StreamingDMD(max_rank) + % StreamingDMD Create an object for incrementaly updated DMD. + % sdmd = StreamingDMD or + % sdmd = StreamingDMD(0) Use the direct algorithm + % (equivalent to batch-processed DMD) + % + % obj = StreamingDMD(max_rank) Use a POD-truncation + % version of the algorithm, with max_rank corresponding + % to the truncation threshold (i.e., the maximum number + % of POD modes retained). + + if nargin == 1 + obj.max_rank = max_rank; + end + end + + function obj = update(obj, x, y) + % UPDATE Update the DMD computation with a pair of snapshots + % + % sdmd = sdmd.UPDATE(x,y) adds the pair of snapshots (x,y) + % to the data ensemble. Here, if the (discrete-time) + % dynamics are given by z(n+1) = f(z(n)), then (x,y) + % should be measurements correponding to consecutive states + % z(n) and z(n+1). + + % parameters + ngram = 5; % number of times to reapply Gram-Schmidt + epsilon = eps; % tolerance for expanding the bases + + obj.count = obj.count + 1; + + % compute snapshot norms + normx = norm(x); + normy = norm(y); + + % ---- Process the First Iterate ---- + if obj.count == 1 + % construct bases + obj.Qx = x / norm(x); + obj.Qy = y / norm(y); + + % compute + obj.Gx = normx^2; + obj.Gy = normy^2; + obj.A = normx * normy; + return + end + + % ---- ALGORITHM STEP 1 ---- + % classical Gram-Schmidt reorthonormalization + xtilde = zeros(size(obj.Qx,2),1); + ytilde = zeros(size(obj.Qy,2),1); + ex = x; + ey = y; + for igram = 1:ngram + dx = obj.Qx' * ex; + dy = obj.Qy' * ey; + xtilde = xtilde + dx; + ytilde = ytilde + dy; + ex = ex - obj.Qx * dx; + ey = ey - obj.Qy * dy; + end + + % ---- ALGORITHM STEP 2 ---- + % check basis for x and expand, if necessary + if norm(ex) / normx > epsilon + % update basis for x + obj.Qx = [obj.Qx ex/norm(ex)]; + % increase size of Gx and A (by zero-padding) + obj.Gx = [obj.Gx zeros(size(obj.Gx,1),1); zeros(1,size(obj.Gx,2)+1)]; + obj.A = [obj.A zeros(size(obj.A,1),1)]; + end + + % check basis for y and expand, if necessary + if norm(ey) / normy > epsilon + % update basis for y + obj.Qy = [obj.Qy ey/norm(ey)]; + % increase size of Gy and A (by zero-padding) + obj.Gy = [obj.Gy zeros(size(obj.Gy,1),1); zeros(1,size(obj.Gy,2)+1)]; + obj.A = [obj.A; zeros(1,size(obj.A,2))]; + end + + % ---- ALGORITHM STEP 3 ---- + % check if POD compression is needed + r0 = obj.max_rank; + if r0 + if size(obj.Qx,2) > r0 + [evec, eval] = eig(obj.Gx); + [eval, indx] = sort(diag(eval),'descend'); + qx = evec(:,indx(1:r0)); + obj.Qx = obj.Qx * qx; + obj.A = obj.A * qx; + obj.Gx = diag(eval(1:r0)); + end + if size(obj.Qy,2) > r0 + [evec, eval] = eig(obj.Gy); + [eval, indx] = sort(diag(eval),'descend'); + qy = evec(:,indx(1:r0)); + obj.Qy = obj.Qy * qy; + obj.A = qy' * obj.A; + obj.Gy = diag(eval(1:r0)); + end + end + + % ---- ALGORITHM STEP 4 ---- + xtilde = obj.Qx' * x; + ytilde = obj.Qy' * y; + + % update A and Gx + obj.A = obj.A + ytilde * xtilde'; + obj.Gx = obj.Gx + xtilde * xtilde'; + obj.Gy = obj.Gy + ytilde * ytilde'; + end + + function [modes, evals] = compute_modes(obj) + % COMPUTE_MODES Compute DMD modes and eigenvalues + % + % [V, lam] = sdmd.COMPUTE_MODES produces a vector lam + % of DMD eigenvalues, and a matrix V whose columns are + % DMD modes of the current dataset. + + Ktilde = obj.Qx' * obj.Qy * obj.A * pinv(obj.Gx); + [evecK, evals] = eig(Ktilde); + evals = diag(evals); + modes = obj.Qx * evecK; + end + end +end \ No newline at end of file diff --git a/2DLinearTimeVarying/WindowDMD.m b/2DLinearTimeVarying/WindowDMD.m new file mode 100644 index 0000000..d82bda1 --- /dev/null +++ b/2DLinearTimeVarying/WindowDMD.m @@ -0,0 +1,155 @@ +% WindowDMD is a class that implements window dynamic mode decomposition +% The time complexity (multiply-add operation for one iteration) is +% O(8n^2), and space complexity is O(2wn+2n^2), where n is the state +% dimension, w is the window size. +% +% Algorithm description: +% At time step k, define two matrix +% X(k) = [x(k-w+1),x(k-w+2),...,x(k)], Y(k) = [y(k-w+1),y(k-w+2), +% ...,y(k)], that contain the recent w snapshot pairs from a finite +% time window, where x(k), y(k) are the n dimensional state vector, +% y(k) = f(x(k)) is the image of x(k), f() is the dynamics. +% Here, if the (discrete-time) dynamics are given by z(k) = f(z(k-1)) +% , then x(k), y(k) should be measurements corresponding to +% consecutive states z(k-1) and z(k). +% At time step k+1, we need to forget old snapshot pair xold = +% x(k-w+1), yold = y(k-w+1), and remember new snapshot pair xnew = +% x(k+1), ynew = y(k+1). +% We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively +% by efficient rank-2 updating window DMD algorithm. +% An exponential weighting factor can be used to place more weight on +% recent data. +% +% Usage: +% wdmd = WindowDMD(n,w,weighting) +% wdmd.initialize(Xw,Yw) +% wdmd.update(xnew, ynew) +% [evals, modes] = wdmd.computemodes() +% +% properties: +% n: state dimension +% w: finite time window size +% weighting: weighting factor in (0,1] +% timestep: number of snapshot pairs processed +% Xw: recent w snapshots x stored in Xw, size n by w +% Yw: recent w snapshots y stored in Yw, size n by w +% A: DMD matrix for w snapshot pairs, size n by n +% P: Matrix that contains information about recent w snapshots, +% size n by n +% +% methods: +% initialize(Xw, Yw), initialize window DMD algorithm +% update(xnew, ynew), update by forgetting old snapshot pairs, +% and remeber new snapshot pair, move sliding window forward +% At time k+1, X(k+1) = [x(k-w+2),x(k-w+2),...,x(k+1)], +% Y(k+1) = [y(k-w+2),y(k-w+2),...,y(k+1)], +% we should take xnew = x(k+1), ynew = y(k+1) +% computemodes(), compute and return DMD eigenvalues and DMD modes +% +% Authors: +% Hao Zhang +% Clarence W. Rowley +% +% Reference: +% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, +% ``Online Dynamic Mode Decomposition for Time-varying Systems", +% in production, 2017. Available on arXiv. +% +% Created: +% April 2017. +% +% To look up the documentation, type help WindowDMD + +classdef WindowDMD < handle + properties + n = 0; % state dimension + w = 0; % window size + weighting = 1; % weighting factor + timestep = 0; % number of snapshots processed + Xw; % recent w snapshots x stored in matrix Xw + Yw; % recent w snapshots y stored in matrix Yw + A; % DMD matrix for w snapshot pairs, size n by n + P; % Matrix that contains information about recent w snapshots + % , size n by n + end + + methods + function obj = WindowDMD(n, w, weighting) + % Creat an object for window DMD + % Usage: wdmd = WindowDMD(n,w,weighting) + if nargin == 3 + obj.n = n; + obj.w = w; + obj.weighting = weighting; + obj.Xw = zeros(n,w); + obj.Yw = zeros(n,w); + obj.A = zeros(n,n); + obj.P = zeros(n,n); + end + end + + function initialize(obj, Xw, Yw) + % Initialize WindowDMD with w snapshot pairs stored in (Xw, Yw) + % Usage: wdmd.initialize(Xw,Yw) + + % initialize Xw, Yw + obj.Xw = Xw; obj.Yw = Yw; + % initialize A, P + q = length(Xw(1,:)); + if(obj.timestep == 0 && obj.w == q && rank(Xw) == obj.n) + weight = (sqrt(obj.weighting)).^(q-1:-1:0); + Xw = Xw.*weight; + Yw = Yw.*weight; + obj.A = Yw*pinv(Xw); + obj.P = inv(Xw*Xw')/obj.weighting; + end + obj.timestep = obj.timestep + q; + end + + function update(obj, xnew, ynew) + % Update the DMD computation by sliding the finite time window + % forward. + % Forget the oldest pair of snapshots (xold, yold), and + % remembers the newest pair of snapshots (xnew, ynew) in the + % new time window. If the new finite time window at time step + % k+1 includes recent w snapshot pairs as + % X(k+1) = [x(k-w+2),x(k-w+3),...,x(k+1)], + % Y(k+1) = [y(k-w+2),y(k-w+3),...,y(k+1)], + % where y(k) = f(x(k)) and f is the dynamics, then we should + % take xnew = x(k+1), ynew = y(k+1) + % Usage: wdmd.update(xnew, ynew) + + % define old snapshots to be discarded + xold = obj.Xw(:,1); yold = obj.Yw(:,1); + % Update recent w snapshots + obj.Xw = [obj.Xw(:,2:end), xnew]; + obj.Yw = [obj.Yw(:,2:end), ynew]; + + % direct rank-2 update + % define matrices + U = [xold, xnew]; V = [yold, ynew]; + C = diag([-(obj.weighting)^(obj.w),1]); + % compute PkU matrix vector product beforehand + PkU = obj.P*U; + % compute AkU matrix vector product beforehand + AkU = obj.A*U; + % compute Gamma + Gamma = inv(inv(C)+U'*PkU); + % update A + obj.A = obj.A + (V-AkU)*(Gamma*PkU'); + % update P + obj.P = (obj.P - PkU*(Gamma*PkU'))/obj.weighting; + % ensure P is SPD by taking its symmetric part + obj.P = (obj.P+(obj.P)')/2; + + % time step + 1 + obj.timestep = obj.timestep + 1; + end + + function [evals, modes] = computemodes(obj) + % Compute DMD eigenvalues and DMD modes at current time step + % Usage: [evals, modes] = wdmd.computemodes() + [modes, evals] = eig(obj.A, 'vector'); + end + end +end \ No newline at end of file