Commit 25807034 authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

new: eig for quantity.Discrete + unittests. implemented pointwise in space /...

new: eig for quantity.Discrete + unittests. implemented pointwise in space / on domain. only for ndims(domain)==1
parent e37e9a54
......@@ -1834,6 +1834,43 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
'name', "exp(" + obj(1).name + ")");
end % exp()
function [V, D, W] = eig(obj)
% eig Eigenvalues and eigenvectors pointwise over the domain
% E = eig(A) produces a column vector E containing the eigenvalues of
% a square matrix A.
%
% [V,D] = eig(A) produces a diagonal matrix D of eigenvalues and
% a full matrix V whose columns are the corresponding eigenvectors
% so that A*V = V*D.
%
% [V,D,W] = eig(A) also produces a full matrix W whose columns are the
% corresponding left eigenvectors so that W'*A = D*W'.
assert(numel(obj(1).domain)==1, "quantity.Discrete/eig is only implemented for one domain");
assert((ndims(obj) <= 2) && (size(obj, 1) == size(obj, 2)), ...
"quantity.Discrete/eig is only implemented for quadratic quantities");
objDisc = permute(obj.on(), [2, 3, 1]);
if nargout == 1
% E = eig(A)
VDisc = zeros([obj(1).domain.n, size(obj, 1)]);
for it = 1 : obj(1).domain.n
VDisc(it, :) = eig(objDisc(:, :, it));
end % for it = 1 : obj(1).domain.n
else
% [V,D] = eig(A) or [V,D,W] = eig(A)
VDisc = zeros([obj(1).domain.n, size(obj)]);
DDisc = zeros([obj(1).domain.n, size(obj)]);
WDisc = zeros([obj(1).domain.n, size(obj)]);
for it = 1 : obj(1).domain.n
[VDisc(it, :, :), DDisc(it, :, :), WDisc(it, :, :)] = ...
eig(objDisc(:, :, it));
end % for it = 1 : obj(1).domain.n
D = quantity.Discrete(DDisc, obj(1).domain, "name", "D");
W = quantity.Discrete(WDisc, obj(1).domain, "name", "W");
end
V = quantity.Discrete(VDisc, obj(1).domain, "name", "V");
end % eig()
function y = log(obj)
% log Natural logarithm.
% log(X) is the natural logarithm of the elements of X.
......@@ -1964,7 +2001,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
p = a(1).gridSize();
q = p(2);
p = p(1);
A = a.on();
obj = a.on();
B = b.on();
% dimensions
......@@ -1977,7 +2014,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
for k = 1 : n % rows of P
for l = 1 : m % columns of P
for r = 1 : o % rows of B or columns of A
P(:, :, k, l) = P( :, :, k, l ) + A( :, :, k, r) .* B( :, :, r, l );
P(:, :, k, l) = P( :, :, k, l ) + obj( :, :, k, r) .* B( :, :, r, l );
end
end
end
......@@ -2068,7 +2105,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
q = p(2);
p = p(1);
A = a.on();
obj = a.on();
B = repmat(b.on(), 1, 1, 1, q);
B = permute(B, [1, 4, 2, 3]);
......@@ -2082,7 +2119,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
for k = 1 : n % rows of P
for l = 1 : m % columns of P
for r = 1 : o % rows of B or columns of A
P(:, :, k, l) = P( :, :, k, l ) + A( :, :, k, r) .* B( :, :, r, l );
P(:, :, k, l) = P( :, :, k, l ) + obj( :, :, k, r) .* B( :, :, r, l );
end
end
end
......
......@@ -4,6 +4,33 @@ function [tests] = testDiscrete()
tests = functiontests(localfunctions);
end
function testEigScalar(tc)
z = quantity.Domain("z", linspace(-1, 1, 11));
[Vz, Dz, Wz] = z.Discrete.eig();
tc.verifyEqual(z.grid, Dz.on());
tc.verifyEqual(ones(z.n, 1), Vz.on());
tc.verifyEqual(ones(z.n, 1), Wz.on());
end % testEigScalar()
function testEigMat(tc)
z = quantity.Domain("z", linspace(-1, 1, 11));
zSym = sym("z");
quan = quantity.Discrete(quantity.Symbolic([1+zSym^2, 2; -zSym, zSym], z));
% all output arguments
[Vq, Dq, Wq] = quan.eig();
tc.verifyEqual(diag(eig(quan.at(0.4))), Dq.at(0.4), "AbsTol", 1e-15);
tc.verifyEqual(on(quan * Vq), on(Vq * Dq), "AbsTol", 1e-15);
tc.verifyEqual(on(Wq' * quan), on(Dq * Wq'), "AbsTol", 1e-15);
% 2 output arguments
[Vq, Dq] = quan.eig();
tc.verifyEqual(diag(eig(quan.at(0.4))), Dq.at(0.4), "AbsTol", 1e-15);
tc.verifyEqual(on(quan * Vq), on(Vq * Dq), "AbsTol", 1e-15);
% 1 output argument
Eq = quan.eig();
tc.verifyEqual(eig(quan.at(0.4)), Eq.at(0.4), "AbsTol", 1e-15);
tc.verifyEqual(on(Eq), on(Dq.diag2vec()), "AbsTol", 1e-15);
end % testEigMat()
function testPower(tc)
%% scalar case
z = quantity.Domain("z", linspace(0, 1, 11));
......
Supports Markdown
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