Added TimeDelayOperator to compute functions with space dependent time delays f(t + g(z))

+quantity/TimeDelayOperator.m

+classdef TimeDelayOperator < handle & matlab.mixin.Copyable
+	%TIMEDELAYOPERATOR class to describe an operator of the form
+	%	D[h](z,t) = h( t + coefficient(z) )
+	% This comes from the inverse Laplace transfromed value of
+	%	exp(-s coefficient(z) ) * H(s)   <->   h(t + coefficient(z))
+	%
+	
+	properties
+		coefficient {mustBe.quantity};
+	end
+	
+	properties (Access = protected)
+		spatialDomain_inner quantity.Domain;
+		timeDomain_inner quantity.Domain;
+		isZero logical = false;
+	end
+	
+	methods
+		function obj = TimeDelayOperator(coefficient, timeDomain, optionalArgs)
+			%TIMEDELAYOPERATOR Construct an instance of a time delay
+			%operator
+			%   obj = TimeDelayOperator(COEFFICIENT, TIMEDOMAIN, VARARGIN)
+			%   will construct a time delay operator of the form
+			%		OBJ[h](z,t) = h( t + COEFFICIENT(z) )
+			%	with the time delay described by COEFFICIENT on a temporal
+			%	domain TIMEDOMAIN. Using the name-value-pairs optional
+			%	parameters, a spatial domain can be specified by the name
+			%	'spatialDomain'.
+			arguments
+				coefficient;
+				timeDomain;
+				optionalArgs.spatialDomain =  coefficient.domain;
+				optionalArgs.isZero logical = false;
+			end
+			
+			obj.coefficient = coefficient;
+			obj.timeDomain_inner = timeDomain;
+			obj.spatialDomain_inner = optionalArgs.spatialDomain;
+			obj.isZero = optionalArgs.isZero;
+		end
+		
+		function d = applyTo(obj, h)
+			%APPLYON apply the operator on a function
+			%   D = applyOn(OBJ, H) computes the application of the
+			%   operator OBJ to the function H. 
+			%		d = OBJ[h]
+			arguments
+				obj
+				h quantity.Discrete;
+			end
+			
+			n = size(obj, 1);
+			m = size(obj, 2); assert(m == size(h, 1));
+			l = size(h, 2);
+			
+			d = quantity.Discrete.zeros([n, l], [obj.spatialDomain, obj.timeDomain]);
+			
+			for i = 1 : n
+				for j = 1 : m
+					for k = 1 : l
+						if ~obj(i, j).isZero
+							d(i, k) = d(i, k) + h(j, k).compose( obj.t + obj(i, j).coefficient );
+						end
+					end
+				end
+			end
+		end		
+		
+		function d = compositionDomain( obj, h)
+			
+			% todo assert all coefficients have the same domain
+			
+			d = h(1).compositionDomain( obj.t + obj(1).coefficient );
+		end
+		
+		function t = t(obj)
+			t = quantity.Symbolic(sym(obj.timeDomain.name), 'domain', obj.timeDomain);
+		end
+		
+		function d = diag(obj)
+			
+			n = length(obj);
+			idx = 1:n;
+			
+			d(1:(n^2)) = obj.zero(obj.timeDomain);
+			d = reshape(d, n, n);
+			d(idx == idx') = obj(:);
+						
+		end
+		
+		function C = coefficients(obj)		
+			C = reshape([obj.coefficient], size(obj));
+		end
+		
+		function sd = spatialDomain(obj)
+			sd = join( [obj.spatialDomain_inner] );
+		end
+		
+		function td = timeDomain(obj)
+			td = join( [obj.timeDomain_inner] );
+		end
+		
+
+	end
+	
+	methods (Static)
+		function o = zero(timeDomain, namedArgs)
+			
+			arguments
+				timeDomain quantity.Domain;
+				namedArgs.spatialDomain quantity.Domain = quantity.Domain.empty();
+			end
+			
+			if isempty( namedArgs.spatialDomain )			
+				infValue = inf;
+			else
+				infValue = inf(obj.spatialDomain.gridLength, 1);
+			end
+			
+			qinf = quantity.Discrete(infValue, ...
+				'domain', namedArgs.spatialDomain);
+			o = quantity.TimeDelayOperator(qinf, timeDomain, 'isZero', true);			
+		end		
+	end
+	
+end
+

+unittests/+quantity/testDiscrete.m

@@ -1093,10 +1093,7 @@ testCase.verifyEqual( TZX1.on(), TZ1.on(), 'AbsTol', 1e-12 );
 testCase.verifyEqual( on( T + 0.5 + X  ), on( subs( TZX, 'z', 0.5) ), 'AbsTol', 1e-12 );
 
 XTZ = X+T+Z;
-
-TZX + XTZ;
-
-%TODO
+testCase.verifyEqual( on( TZX + XTZ ), on( 2 * TZX ), 'AbsTol', eps )
 
 end
 

+unittests/+quantity/testTimeDelayOperator.m

+function [tests] = testTimeDelayOperator()
+%TESTGEVREY Summary of this function goes here
+%   Detailed explanation goes here
+tests = functiontests(localfunctions);
+end
+
+function setupOnce(testCase)
+
+z = quantity.Domain('grid', linspace(0,1), 'name', 'z');
+zeta = quantity.Domain('grid', linspace(0,1), 'name', 'zeta');
+t = quantity.Domain('grid', linspace(0, 2, 31), 'name', 't');
+v = quantity.Discrete( z.grid, 'domain', z );
+c_zeta = quantity.Discrete( zeta.grid * 0, 'domain', zeta );
+c = quantity.Discrete( 1, 'domain', quantity.Domain.empty);
+
+testCase.TestData.z = z;
+testCase.TestData.t = t;
+testCase.TestData.zeta = zeta;
+testCase.TestData.coefficient.fun_zeta = c_zeta;
+testCase.TestData.coefficient.variable = v;
+testCase.TestData.coefficient.constant = c;
+	
+testCase.TestData.delay.variable = quantity.TimeDelayOperator(v, t);
+testCase.TestData.delay.constant = quantity.TimeDelayOperator(c, t);
+testCase.TestData.delay.zero = quantity.TimeDelayOperator.zero(t);
+testCase.TestData.delay.spatialDomain2 = quantity.TimeDelayOperator( ...
+	v - c_zeta, t);
+
+testCase.TestData.delay.diagonal = diag( ...
+	[testCase.TestData.delay.constant, ...
+	 testCase.TestData.delay.variable]);
+end
+
+function testInit(testCase)
+
+td = testCase.TestData;
+
+testCase.verifyEqual(td.delay.constant.coefficient, ...
+	td.coefficient.constant);
+testCase.verifyEqual( td.delay.variable.coefficient, td.coefficient.variable);
+	
+end
+
+
+function testDiag(testCase)
+
+	D = testCase.TestData.delay.diagonal;
+	v = testCase.TestData.delay.variable;
+	c = testCase.TestData.delay.constant;
+	o = testCase.TestData.delay.zero;
+
+	testCase.verifyEqual( D(1), c );
+	testCase.verifyEqual( D(2), o );
+	testCase.verifyEqual( D(3), o );
+	testCase.verifyEqual( D(4), v );
+	
+end
+
+function testApplyTo(testCase)
+
+	z = testCase.TestData.z;
+	t = testCase.TestData.t;
+	tau = quantity.Domain('grid', linspace(-1, 3, 201), 'name', 'tau');
+	v = testCase.TestData.delay.variable;
+	
+	h = quantity.Discrete( sin(tau.grid * 2 * pi), 'domain', tau );
+	H = quantity.Discrete( sin(z.grid * 2 * pi + t.grid' * 2 * pi), ...
+		'domain', [z, t]);
+	H_ = v.applyTo(h);
+	
+	%% test a scalar operator
+	testCase.verifyEqual(H.on(), H_.on(), 'AbsTol', 5e-3);
+	
+	d2 = testCase.TestData.delay.spatialDomain2;
+	H2 = d2.applyTo(h);
+	H2 = subs(H2, 'zeta', 0);
+	testCase.verifyEqual(H.on(), H2.on(), 'AbsTol', 5e-3);
+	
+	% test a diagonal operator
+	h = quantity.Discrete( {sin(tau.grid * 2 * pi); ...
+		sin(tau.grid * 2 * pi)}, 'domain', tau );
+	H = quantity.Discrete( {sin(z.grid * 0 + 2 * pi + t.grid' * 2 * pi); ...
+		sin(z.grid * 2 * pi + t.grid' * 2 * pi)}, 'domain', [z, t]);
+	D = testCase.TestData.delay.diagonal;
+	H_ = D.applyTo(h);
+	
+	testCase.verifyEqual( H.on(), H_.on(), 'AbsTol', 5e-3)
+	
+	
+
+	
+end
+
+
+
+
+
+
+
+
+
+