Minor changes and bug fixes

+misc/formatMatrix.m

+function [s] = formatMatrix(M)
+
+s = '';
+
+for k = 1:size(M,1)
+	s = [s, sprintf(' %.3g', M(k,:)), '\n'];
+end
+
+end
+

+numeric/near.m

-function [i, m] = near(nominalValue, actualValue, tolerance)
+function [i, m, s] = near(nominalValue, actualValue, tolerance, relative)
 %NEAR checks 'equality' for double numerical values.
 %   i = near(nominalValue, actualValue, tolerance) checks if.
 %       abs(nominalValue - actualValue) < tolerance.
@@ -12,12 +12,23 @@ if nargin == 1
 	actualValue = zeros(size(nominalValue));
 end
 
-absV = abs(nominalValue - actualValue);
+if nargin < 4
+	relative = 1;
+else
+	if islogical(relative)
+		relative = max(abs(nominalValue(:)));
+	end
+end
+	
+
+absV = abs(nominalValue ./ relative - actualValue ./ relative);
 i = all(absV(:) < tolerance);
 m = max( absV(:) );
 
+s = sprintf('Is near %i, maximal difference: %d', [i, m]);
+
 if nargout == 0
-	i = sprintf('Is near %i, maximal difference: %d', [i, m]);
+	i = s;
 end
 
 end
\ No newline at end of file

+quantity/Operator.m

@@ -149,10 +149,24 @@ classdef Operator < handle & matlab.mixin.Copyable
 		function C = uminus(A)
 			C = (-1)*A;
 		end
-		
+		function s = size(obj, varargin)
+			s = size(obj(1).coefficient, varargin{:});
+		end
+		function l = length(obj)
+			l = numel(obj);
+		end
 	end
 	
 	methods (Access = public)
+		function d = double(obj)
+			d = double(obj.M);
+		end
+		function [i, m, s] = near(obj, B, varargin)
+			[i, m, s] = obj.M.near(B.M, varargin{:});
+			if nargout == 0
+				i = s;
+			end
+		end
 		function i = isempty(obj)
 			
 			i = numel(obj)==0 || isempty(obj(1).coefficient);
@@ -240,8 +254,8 @@ classdef Operator < handle & matlab.mixin.Copyable
 			spatialDomain = A(1).coefficient(1).grid{:};
 			
 			% build the large system of odes from the recursion
-			A_full = zeros(length(spatialDomain), d*n, d*n);
-			B_full = zeros(length(spatialDomain), d*n, m);
+			A_full = zeros(length(spatialDomain), N*n, N*n);
+			B_full = zeros(length(spatialDomain), N*n, m);
 			
 			idx = 1:n;
 			for k = 1:N
@@ -249,9 +263,9 @@ classdef Operator < handle & matlab.mixin.Copyable
 				if length(B) >= k
 					B_full(:, idxk, :) = B(k).coefficient.on();
 				end
-				for j = 1:d
+				for j = 1:N
 					idxj = idx + (j-1)*n;
-					if k - (j-1) > 0
+					if k - (j-1) > 0 && k - (j-1) <= length(A)
 						A_full(:, idxk, idxj) = A(k - (j-1)).coefficient.on();
 					end
 				end
@@ -261,11 +275,12 @@ classdef Operator < handle & matlab.mixin.Copyable
 				'grid', {spatialDomain, spatialDomain}, 'gridName', {'z', 'zeta'});
 			
 			p = f.volterra(quantity.Discrete(B_full, 'grid', spatialDomain, ...
-				'gridName', 'zeta'));
+				'gridName', 'zeta'), 'variable', 'zeta');
 			
+			z0 = A(1).grid{1}(1);
 			for k = 1:N
 				idxk = idx + (k-1)*n;
-				Phi(:,:,k) = f(idxk, idx).subs('zeta', 0);
+				Phi(:,:,k) = f(idxk, idx).subs('zeta', z0);
 				Psi(:,:,k) = p(idxk,:);
 			end
 			
@@ -331,7 +346,7 @@ classdef Operator < handle & matlab.mixin.Copyable
 			%   Psi_{k-i}(z, xi) + B_k(xi) d_xi
 			
 			
-			F0 = A.computePhi0();
+			F0 = quantity.Discrete( A.computePhi0() );
 			Phi(:,:,1) = F0.subs('zeta', 0);
 			Psi(:,:,1) = F0.volterra(B(1).coefficient.subs('z', 'zeta'));
 			
@@ -388,7 +403,7 @@ classdef Operator < handle & matlab.mixin.Copyable
 					% exponential function:
 					z = sym('z', 'real');
 					zeta = sym('zeta', 'real');
-					f0 = expm(squeeze(obj(1).coefficient.on())*(z - zeta));
+					f0 = expm(obj(1).coefficient.atIndex(1)*(z - zeta));
 					F0 = quantity.Symbolic(f0, 'grid', {spatialDomain, spatialDomain});
 				else
 					f0 = misc.fundamentalMatrix.odeSolver_par( ...
@@ -409,6 +424,9 @@ classdef Operator < handle & matlab.mixin.Copyable
 		function newOperator = changeGrid(obj, gridNew, gridNameNew)
 			newOperator = quantity.Operator(obj.M.changeGrid(gridNew, gridNameNew));
 		end
+		function newOperator = subs(obj, varargin)
+			newOperator = obj.M.subs(varargin{:});
+		end
 	end
 	
 	methods (Static, Access = protected)

+quantity/Settings.m

 classdef Settings < handle
 	%SETTINGS object to save some settings for the quantity class
 	%   
-	
+
 	properties (Access = public)
-		dockThePlot = false
+		plot = struct('dock', true, 'showTitle', true);
+		progressBar = struct('show', true);
 	end
 	
 	methods (Access = private)

+unittests/+quantity/testDiscrete.m

@@ -996,39 +996,6 @@ testCase.verifyEqual(referenceResult, permute(ax2.on(), [1, 3, 2]), 'AbsTol', 1e
 %%
 end
 
-function testMTimes2(testCase)
-
-%% TODO@ff
-% z = sym('z', 'real');
-% zeta = sym('zeta', 'real');
-%
-% Z = linspace(0,1)';
-% A = [ 0, 0, 0, -3; ...
-% 	0, 0, 0,  0; ...
-% 	0, 0, 0,  0; ...
-% 	0, 0, 0,  0];
-%
-% B = [              1,              0,        0, 0; ...
-% 	z - zeta,              1,        0, 0; ...
-% 	(z - zeta)^2/2,       z - zeta,        1, 0; ...
-% 	(z - zeta)^3/6, (z - zeta)^2/2, z - zeta, 1];
-%
-% a = quantity.Symbolic(sym(A), 'grid', {Z, Z'}, 'variable', [z, zeta]);
-% b = quantity.Symbolic(B, 'grid', {Z, Z'}, 'variable', [z, zeta]);
-%
-% ab = A * B;
-%
-% p = quantity.Discrete(a);
-% q = quantity.Discrete(b);
-%
-% pq = matTimes( p, q );
-% [pq.figureID] = deal(2);
-%
-% %
-% verifyTrue(testCase, numeric.near(pq.on(), ab.on));
-
-end
-
 function testZeros(testCase)
 z = quantity.Discrete.zeros([2, 7, 8], {linspace(0,10)', linspace(0, 23, 11)}, 'gridName', {'a', 'b'});
 O = zeros(100, 11, 2, 7, 8);

+unittests/+quantity/testFunction.m

@@ -4,6 +4,17 @@ function [tests] = testFunction()
 tests = functiontests(localfunctions);
 end
 
+function testOn(testCase)
+z = linspace(0, 2*pi)';
+f = quantity.Function({@(z) sin(z), @(z) cos(z)}, 'grid', z, 'gridName', 'z');
+
+F1 = f.on(z);
+F2 = f.on();
+
+testCase.verifyTrue( numeric.near(F1, F2) );
+
+end
+
 function testUMinus(testCase)
 f = quantity.Function(@(z) z, 'grid', linspace(0,1).');
 mf = -f;
@@ -66,9 +77,9 @@ f = quantity.Function({f1 ; f2 }, 'grid', linspace(0,1,1e4).');
 
 F = f.inner(f.');
 
-verifyEqual(testCase, integral(@(z) f1(z).^2, 0, 1), F(1,1));
-verifyEqual(testCase, integral(@(z) f2(z).^2, 0, 1 ), F(2, 2));
-verifyEqual(testCase, integral(@(z) f1(z) .* f2(z), 0, 1 ) , F(1, 2));
+verifyTrue(testCase, numeric.near( integral(@(z) f1(z).^2, 0, 1), F(1,1), 1e-6));
+verifyTrue(testCase, numeric.near( integral(@(z) f2(z).^2, 0, 1 ), F(2, 2), 1e-6));
+verifyTrue(testCase, numeric.near( integral(@(z) f1(z) .* f2(z), 0, 1 ) , F(1, 2), 1e-6));
 end
 
 function testMTimes(testCase)
@@ -115,21 +126,6 @@ verifyEqual(testCase, f.gridSize, 42);
 
 end
 
-function testPlot(testCase)
-%%
-f = quantity.Function(...
-    {@(z) (sin(z .* pi) ), @(z) z; @(z) 0.*z, @(z)(cos(z .* pi) )}, ...
-    'grid', {linspace(0, 1, 42).'}, ...
-    'gridName', {'z'});
-
-fig = f.plot('fig', 23);
-
-%%
-verifySameHandle(testCase, fig, gcf);
-close(fig);
-end
-
-
 function testEvaluate(testCase)
 %%
 f1 = @(z) (sin(z .* pi) );

+unittests/+quantity/testOperator.m

@@ -37,20 +37,6 @@ testCase.TestData.a = a;
 
 end
 
-function testPlot(testCase)
-	f = testCase.TestData.A.plot();
-	close(f);
-end
-
-function testPowerSeries(testCase)
-	syms z s
-	P = testCase.TestData.A.powerSeries();
-	testCase.verifyEqual(P(1, 4), - 4*s^2 + (z - 3)*s);
-	
-	error('Write test for numerical Values')
-	
-end
-
 function testCTranspose(testCase)
 	At = testCase.TestData.A';
 	testCase.verifyTrue(At(3).coefficient == testCase.TestData.A(3).coefficient');	
@@ -77,14 +63,14 @@ function testFundamentalMatrixSpaceDependent(testCase)
 
 	PhiExact = quantity.Symbolic([exp(z), f; 0, exp(a*z)], 'grid', Z, 'name', 'Phi_A');
 		
-	testCase.verifyTrue(numeric.near(F0.on() ./ F0.MAX(), PhiExact.on() ./ F0.MAX(), 1e-6));
+	testCase.verifyTrue(numeric.near(F0.on(), PhiExact.on(), 1e-6, true));
 	
 end
 % 
 function testStateTransitionMatrix(testCase)
-N = 3;
+N = 4;
 z = sym('z', 'real');
-Z = linspace(0, 1, 51)';
+Z = linspace(0, 1, 11)';
 a = quantity.Symbolic(cat(3, ...
    [ 0, 0, 0, 0; ...
      1   0  0  0; ...
@@ -100,34 +86,26 @@ a = quantity.Symbolic(cat(3, ...
     0 0 0 0]...
    ), 'name', 'A', 'grid', Z);
 A = quantity.Operator(a);
-B = quantity.Symbolic([1 1; zeros(3,2)], 'name', 'B', 'variable', z, 'grid', Z);
+B = quantity.Symbolic([1 0; 0 0; 0 1; 0 0], 'name', 'B', 'variable', z, 'grid', Z);
 B = quantity.Operator(B);
 
 [Phi1, Psi1] = A.stateTransitionMatrix(N, B);
 [Phi2, Psi2] = A.stateTransitionByOde(N, B);
 
-testCase.verifyTrue(numeric.near(Phi1.M.on(), Phi2.M.on(), 1e-3));
-testCase.verifyTrue(numeric.near(Psi1.M.on(), Psi2.M.on(), 1e-3));
+testCase.verifyTrue(numeric.near(Phi1.M.on(), Phi2.M.on(), 1e-2));
+testCase.verifyTrue(numeric.near(Psi1.M.on(), Psi2.M.on(), 1e-2));
 
-%% TODO: write a test for the fundamental matrices.
-% gamma = ones(4,1);
-% nu = ones(2,1);
-% 
-% for k = 1:N
-% 	wk_num(k,:) = Phi(:,:,k) * gamma + Psi(:,:,k) * nu;
-% end
-% 	
-% wk_ode = quantity.Discrete.empty(0,4);
-% 
-% for k = 1:N
-% % odeStep = @(z, w, A, B, v) (A.on(z) * w + B.on(z) * v);
-% 	[~, tmp] = ode15s(@(z,w) odeStep(z, w, A, B, nu, wk_ode), Z, gamma * (k==1));
-% 	wk_ode(k, :) = quantity.Discrete(tmp, 'grid', Z, 'gridName', 'z');
-% end
-% 
-% P = wk_num.relativeErrorSupremum(wk_ode);
-% P.plot
-% 
+%% test also the parabolic system
+A0 = quantity.Symbolic([0 1; 1 0], 'grid', Z, 'variable', 'z');
+A1 = quantity.Symbolic([0 1; 0 0], 'grid', Z, 'variable', 'z');
+A = quantity.Operator({A0, A1});
+B = quantity.Operator(quantity.Symbolic([-1 -1; 0 0], 'grid', Z, 'variable', 'z'));
+
+[Phi1, Psi1] = A.stateTransitionMatrix(N, B);
+[Phi2, Psi2] = A.stateTransitionByOde(N, B);
+
+testCase.verifyTrue(Phi1.near(Phi2, 1e-2));
+testCase.verifyTrue(Psi1.near(Psi2, 1e-2));
 
 end
 
@@ -199,7 +177,7 @@ function testMTimes(testCase)
 	Cs = misc.polynomial2coefficients(As*Bs, s);
 	
 	for k = 1:3
-		testCase.verifyTrue( numeric.near(Cs(:,:,k), C(k).coefficient.at(0)) );s
+		testCase.verifyTrue( numeric.near(Cs(:,:,k), C(k).coefficient.at(0)) );
 	end
 	
 	
@@ -208,7 +186,7 @@ function testMTimes(testCase)
 	C =	K*A;
 	
 	for k = 1:3
-		numeric.near(double(C(k).coefficient), double(K * A(k).coefficient))
+		numeric.near(double(C(k).coefficient), double(K * A(k).coefficient));
 	end
 	
 end
@@ -226,16 +204,16 @@ z = linspace(0, pi, 15)';
 	B = quantity.Operator({A0});
 	
 	C = A + B;
-	
-	numeric.near(double(C(1).coefficient), double(A0 + A0))
-	numeric.near(double(C(2).coefficient), double(A1))
-	numeric.near(double(C(3).coefficient), double(A2))
+		
+	testCase.verifyTrue( numeric.near(double(C(1).coefficient), double(A0 + A0)) );
+	testCase.verifyTrue( numeric.near(double(C(2).coefficient), double(A1)) );
+	testCase.verifyTrue( numeric.near(double(C(3).coefficient), double(A2)) );
 	
 	C = B + A;
 	
-	numeric.near(double(C(1).coefficient), double(A0 + A0))
-	numeric.near(double(C(2).coefficient), double(A1))
-	numeric.near(double(C(3).coefficient), double(A2))
+	testCase.verifyTrue( numeric.near(double(C(1).coefficient), double(A0 + A0)) );
+	testCase.verifyTrue( numeric.near(double(C(2).coefficient), double(A1)) );
+	testCase.verifyTrue( numeric.near(double(C(3).coefficient), double(A2)) );
 	
 end