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

quantity.Discrete: simple implementation of det for nxn matrices (shitty for large arrays)

parent 93abf2de
......@@ -2151,8 +2151,12 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
end
end % abs()
function d = det(obj)
function d = det(obj, setName)
% det(X) returns the the determinant of the squre matrix X
arguments
obj;
setName = true;
end
if isempty(obj)
d = quantity.Discrete.empty(1);
elseif numel(obj) == 1
......@@ -2168,12 +2172,21 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
- prod(obj(logical([1, 0, 0; 0, 0, 1; 0, 1, 0]))) ...
- prod(obj(logical([0, 1, 0; 1, 0, 0; 0, 0, 1])));
else
error("only implemented for 2x2 matrices yet");
% maybe a more sophisticated implementation would be nice.
d = 0;
for it = 1 : size(obj, 1)
selector = true(size(obj));
selector(:, 1) = false;
selector(it, :) = false;
d = d + (-1)^(it-1)*obj(it, 1) * det(reshape(obj(selector), size(selector)-1), false);
end % for it = 1 : size(obj, 1)
end
else
error("det is only defined for quadratic matrices");
end
if setName
d.setName("det(" + obj(1).name + ")");
end
end % det()
function y = real(obj)
......
......@@ -15,11 +15,17 @@ end
function testDet(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
zDiscrete = quantity.Discrete(z.grid, 'domain', z);
myMatrix2x2 = [1, 2; 3, 4] * (1 + 0 * zDiscrete);
myMatrix3x3 = [1, 2, 3; 2, 3, 4; 3, 4, 5] * (1 + 0 * zDiscrete);
myMatrix2x2 = [1, 2; 3, 4] * zDiscrete;
myMatrix3x3 = [1, 2, 3; 2, -3, 4; 3, 4, 5] * zDiscrete;
myMatrix5x5 = [1, 2, 3, 4, 0; ...
2, 3, 4, -5, 6; ...
3, 4, 5, 6, 7; ...
4, -6, 7, 8, 9; ...
1, 1, 1, 1, 1] * zDiscrete;
tc.verifyEqual(zDiscrete.det.on(), zDiscrete.on(), 'AbsTol', 10*eps);
tc.verifyEqual(myMatrix2x2.det.at(0.5), det(myMatrix2x2.at(0.5)), 'AbsTol', 10*eps);
tc.verifyEqual(myMatrix3x3.det.at(0.5), det(myMatrix3x3.at(0.5)), 'AbsTol', 10*eps);
tc.verifyEqual(myMatrix5x5.det.at(0.5), det(myMatrix5x5.at(0.5)), 'AbsTol', 10*eps);
end % testDet()
function testSum(tc)
......
......@@ -7,11 +7,17 @@ end
function testDet(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
zSymbolic = quantity.Symbolic(sym("z"), 'domain', z);
myMatrix2x2 = [1, 2; 3, 4] * (1 + 0 * zSymbolic);
myMatrix3x3 = [1, 2, 3; 2, 3, 4; 3, 4, 5] * (1 + 0 * zSymbolic);
myMatrix2x2 = [1, 2; 3, 4] * zSymbolic;
myMatrix3x3 = [1, 2, 3; 2, -3, 4; 3, 4, 5] * zSymbolic;
myMatrix5x5 = [1, 2, 3, 4, 5; ...
2, 3, 4, 5, 6; ...
3, 4, 5, 6, 7; ...
4, -6, 7, 8, 9; ...
1, 1, 1, 1, 1] * zSymbolic;
tc.verifyEqual(zSymbolic.det.on(), zSymbolic.on(), 'AbsTol', 10*eps);
tc.verifyEqual(myMatrix2x2.det.at(0.5), det(myMatrix2x2.at(0.5)), 'AbsTol', 10*eps);
tc.verifyEqual(myMatrix3x3.det.at(0.5), det(myMatrix3x3.at(0.5)), 'AbsTol', 10*eps);
tc.verifyEqual(myMatrix5x5.det.at(0.5), det(myMatrix5x5.at(0.5)), 'AbsTol', 10*eps);
end % testDet()
function testSum(tc)
......
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