Commit 5a460bcd authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

quantity.Discrete: new methods jordanReal & luDecomposition & eye

quantity.Symbolic: new methods ones & eye
parent 25807034
......@@ -1871,6 +1871,53 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
V = quantity.Discrete(VDisc, obj(1).domain, "name", "V");
end % eig()
function [l, u] = luDecomposition(obj)
% luDecomposition split matrix OBJ into lower triangular matrix l and upper triangular
% matrix u
n = size(obj, 1);
assert(numel(obj) == n*n, "obj must be a quadratic matrix");
u = obj;
l = obj.eye([n, n], obj(1).domain);
% simple algorithm from wikipedia:
for it = 1 : n-1
for k = it+1 : n
l(k, it) = u(k, it) / u(it, it);
for jt = it : n
u(k, jt) = u(k, jt) - l(k, it) * u(it, jt);
end % for jt
end % for k
end % for it
end % luDecomposition()
function [V, J, l] = jordanReal(obj)
% jordanReal Real Jordan form of the matrix OBJ pointwise over the (scalar) domain.
% [V, J] = jordanReal(OBJ) computes the transformation matrix V and the jordan
% block matrix J, so that OBJ*V = V*J. The columns of OBJ are the generalised
% eigenvectors.
%
% [V, J, l] = jordanReal(OBJ) also returns the lengths of the jordan
% blocks. l is a vector containing the lengths. So length(l) is the
% number of jordan blocks.
assert(numel(obj(1).domain)==1, "quantity.Discrete/jordanReal is only implemented for one domain");
assert((ndims(obj) <= 2) && (size(obj, 1) == size(obj, 2)), ...
"quantity.Discrete/jordanReal is only implemented for quadratic quantities");
objDisc = permute(obj.on(), [2, 3, 1]);
VDisc = zeros([obj(1).domain.n, size(obj)]);
JDisc = zeros([obj(1).domain.n, size(obj)]);
lDisc = zeros([obj(1).domain.n, size(obj, 1)]);
for it = 1 : obj(1).domain.n
[VDisc(it, :, :), JDisc(it, :, :), lTemp] = ...
misc.jordanReal(objDisc(:, :, it));
lDisc(it, 1:numel(lTemp)) = lTemp;
end % for it = 1 : obj(1).domain.n
V = quantity.Discrete(VDisc, obj(1).domain, "name", "V");
J = quantity.Discrete(JDisc, obj(1).domain, "name", "J");
l = quantity.Discrete(lDisc, obj(1).domain, "name", "l");
end % jordanReal()
function y = log(obj)
% log Natural logarithm.
% log(X) is the natural logarithm of the elements of X.
......@@ -2580,7 +2627,31 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
O = ones([domain.gridLength, valueSize(:)']);
P = quantity.Discrete(O, domain, 'size', valueSize, varargin{:});
end
end
end % ones()
function I = eye(N, domain, varargin)
% eye Identity matrix.
% eye(N, domain) is the N-by-N identity matrix on the domain domain.
%
% eye([M,N], domain) is an M-by-N matrix with 1's on the diagonal and zeros
% elsewhere.
% eye([M,N], domain, varargin) passes further inputs to the quantity.Discrete
% constructor.
if isscalar(N)
N = [N, N];
elseif ~ismatrix(N)
error("eye only supports ndims(N) == 1 or 2");
end
if any( N == 0)
P = quantity.Discrete.empty(N);
else
IvalueDiscrete = reshape(eye(N), [ones(1, numel(domain)), N]);
IvalueDiscrete = repmat(IvalueDiscrete, [domain.n, 1]);
I = quantity.Discrete(IvalueDiscrete, domain, varargin{:});
end
end % eye()
function P = zeros(valueSize, domain, varargin)
%ZEROS initializes an zero quantity.Discrete object
......
......@@ -829,7 +829,25 @@ classdef Symbolic < quantity.Function
%ZEROS initializes an zero quantity.Discrete object of the size
%specified by the input valueSize.
P = quantity.Symbolic(zeros(valueSize), varargin{:});
end
end % zeros()
function P = ones(valueSize, varargin)
%ONES initializes an ones-quantity.Discrete object
% P = ones(VALUESIZE, DOMAIN) creates a matrix of size
% VALUESIZE on the DOMAIN with ones as entries.
P = quantity.Symbolic(ones(valueSize), varargin{:});
end % ones()
function P = eye(N, varargin)
% eye Identity matrix.
% eye(N, domain) is the N-by-N identity matrix on the domain domain.
%
% eye([M,N], domain) is an M-by-N matrix with 1's on the diagonal and zeros
% elsewhere.
% eye([M,N], domain, varargin) passes further inputs to the quantity.Discrete
% constructor.
P = quantity.Symbolic(eye(N), varargin{:});
end % eye()
end % methods (Static)
......
......@@ -4,6 +4,33 @@ function [tests] = testDiscrete()
tests = functiontests(localfunctions);
end
function testEye(tc)
myDomain = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 5))];
tc.verifyEqual(MAX(abs(quantity.Discrete.eye(1, myDomain) - 1)), 0);
tc.verifyEqual(MAX(abs(quantity.Discrete.eye(2, myDomain) - eye(2))), 0);
tc.verifyEqual(MAX(abs(quantity.Discrete.eye([2, 1], myDomain) - eye([2, 1]))), 0);
tc.verifyEqual(MAX(abs(quantity.Discrete.eye([2, 2], myDomain(2)) - eye([2, 2]))), 0);
tc.verifyEqual(MAX(abs(quantity.Discrete.eye([4, 4], myDomain) - eye(4))), 0);
end % testEye()
function testLuDecomposition(tc)
myDomain = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 5))];
syms z zeta
quanSymbolic = quantity.Symbolic([z*zeta, z+1, 2; z^2, zeta^2, -2; 1, 2, 3], myDomain);
quan = quantity.Discrete(quanSymbolic);
% discrete case
[l, u] = quan.luDecomposition();
tc.verifyEqual(MAX(abs(l*u - quan)), 0, 'AbsTol', 1e3*eps);
tc.verifyEqual(MAX(abs(l.*triu(ones(size(l)), 1))), 0, 'AbsTol', 1e3*eps);
tc.verifyEqual(MAX(abs(u.*tril(ones(size(u)), -1))), 0, 'AbsTol', 1e3*eps);
% symbolic case
[l, u] = quanSymbolic.luDecomposition();
tc.verifyEqual(MAX(abs(l*u - quanSymbolic)), 0);
tc.verifyEqual(MAX(abs(l.*triu(ones(size(l)), 1))), 0);
tc.verifyEqual(MAX(abs(u.*tril(ones(size(u)), -1))), 0);
tc.verifyTrue(isa(quanSymbolic, "quantity.Symbolic"));
end % testLuDecomposition()
function testEigScalar(tc)
z = quantity.Domain("z", linspace(-1, 1, 11));
[Vz, Dz, Wz] = z.Discrete.eig();
......@@ -31,6 +58,19 @@ 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 testJordanRealMat(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
[V, J, l] = quan.jordanReal();
[V0p4, J0p4, l0p4] = misc.jordanReal(quan.at(0.4));
tc.verifyEqual(J0p4, J.at(0.4), "AbsTol", 1e-15);
tc.verifyEqual(V0p4, V.at(0.4), "AbsTol", 1e-15);
tc.verifyEqual([l0p4; 0], l.at(0.4), "AbsTol", 1e-15);
tc.verifyEqual(on(quan * V), on(V * J), "AbsTol", 1e-15);
end % testJordanRealMat()
function testPower(tc)
%% scalar case
z = quantity.Domain("z", linspace(0, 1, 11));
......
......@@ -4,6 +4,25 @@ function [tests ] = testSymbolic()
tests = functiontests(localfunctions());
end
function testOnes(tc)
myDomain = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 5))];
tc.verifyEqual(MAX(abs(quantity.Symbolic.ones(1, myDomain) - 1)), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.ones(2, myDomain) - ones(2))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.ones([2, 1], myDomain) - ones([2, 1]))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.ones([2, 2], myDomain(2)) - ones([2, 2]))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.ones([4, 4], myDomain) - ones(4))), 0);
end % testOnes()
function testEye(tc)
myDomain = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 5))];
tc.verifyEqual(MAX(abs(quantity.Symbolic.eye(1, myDomain) - 1)), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.eye(2, myDomain) - eye(2))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.eye([2, 1], myDomain) - eye([2, 1]))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.eye([2, 2], myDomain(2)) - eye([2, 2]))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.eye([4, 4], myDomain) - eye(4))), 0);
end % testEye()
function testValueContinuousBug(tc)
z = quantity.Domain("z", linspace(0, 1, 7));
zeta = quantity.Domain("zeta", linspace(0, 1, 5));
......
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