Skip to content
Snippets Groups Projects
Commit 84ad6074 authored by Felix Kong's avatar Felix Kong
Browse files

compute matrix power Lambda^(k-1) more efficiently

not sure how much my method drifts over the course of thousands of timesteps.
parent 010edfcd
Branches
No related merge requests found
...@@ -94,11 +94,14 @@ classdef dmd < handle ...@@ -94,11 +94,14 @@ classdef dmd < handle
xhat(:,1) = initial_condition; xhat(:,1) = initial_condition;
b0 = obj.Phi\initial_condition; b0 = obj.Phi\initial_condition;
LambdaPowKMinus1 = eye(size(obj.Lambda));
for k = 2:num_timesteps for k = 2:num_timesteps
xhat(:,k) = obj.Phi*obj.Lambda^(k-1)*b0; LambdaPowKMinus1 = obj.Lambda*LambdaPowKMinus1;
end % check = obj.Lambda^(k-1);
% checknorm = norm(LambdaPowKMinus1 - check)/norm(check);
assert(norm(imag(xhat)) <= 1e-6); xhat(:,k) = obj.Phi*LambdaPowKMinus1*b0;
end
% assert(norm(imag(xhat)) <= 1e-6); % assert(norm(imag(xhat)) <= 1e-6);
if (norm(imag(xhat))/norm(xhat) >= 1e-6) % rel norm of imag components if (norm(imag(xhat))/norm(xhat) >= 1e-6) % rel norm of imag components
warning('rel norm of imaginary component of X2hat nontrivial: %e. Proceeding anyway.',norm(imag(xhat))/norm(xhat)) warning('rel norm of imaginary component of X2hat nontrivial: %e. Proceeding anyway.',norm(imag(xhat))/norm(xhat))
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment