made signals.BasicVariable.diff conistent with quantity.Discrete.diff and some...

made signals.BasicVariable.diff conistent with quantity.Discrete.diff and some further beautifications

+fractional/fDiff.m → +fractional/derivative.m

-function [DaY, t] = fDiff(y, p, t0, t, optArg)
-% FDIFF(Y, A, t0, t) calculates the Caputo fractional derivative of
-% order A of a symbolic function Y.
-% ??? assuming that the input was 0 for t <= t0.
-% The output DAY is symbolic.
+function [h, t] = derivative(y, p, t0, t, optArg)
+% DERIVATIVE calculates a fractional derivative
+%	[h, t] = derivative(y, p, t0, t, optArg) comput the p-th fractional derivative of a symbolic
+%	expression y. The lower terminal is t0 and the symbolic independent variable is t. By the
+%	name-value-pair argument, the "type" of the fractional derivative can be specified. Use "RL" for
+%	the Riemann-Liouville derivative and "C" for the Caputo derivative.
 arguments
 	y sym
 	p (1,1) double {mustBePositive};
@@ -16,7 +17,7 @@ if p < 0
 end
 
 if misc.nearInt(p)
-	DaY = diff(y, round(p));
+	h = diff(y, round(p));
 	return
 end
 
@@ -24,12 +25,12 @@ n = ceil(p);
 tau = sym( string(t) + "_");
 
 if optArg.type == "C"	
-	DaY = 1/gamma(n-p) * int((t-tau)^(n-p-1) * subs(diff(y, t, n), t, tau), tau, t0, t);
+	h = 1/gamma(n-p) * int((t-tau)^(n-p-1) * subs(diff(y, t, n), t, tau), tau, t0, t);
 	%DaY = 1/gamma(1+n-a) * ((t - t0)^(n-a) * subs( diff(y, n), t, t0) ...
 	%	+ int((t-tau)^(n-a) * subs(diff(y, 1+n), t, tau), tau, t0, t)); % Integration by parts
 	
 elseif optArg.type == "RL"
-	DaY = 1 / gamma(n-p) * diff( int( (t-tau)^(n-p-1) * subs(y, t, tau), tau, t0, t), t, n);
+	h = 1 / gamma(n-p) * diff( int( (t-tau)^(n-p-1) * subs(y, t, tau), tau, t0, t), t, n);
 else
 	error("The derivative of the type " + optArg.type + " is not implemented. Only C for Caputo and RL for Riemann-Liouville are available")
 end

+signals/BasicVariable.m

@@ -25,13 +25,13 @@ classdef BasicVariable < handle
 			h = numel(obj(1).derivatives);
 		end
 		function T = get.T(obj)
-			T = obj.timeDomain.upper();
+			T = obj.domain.upper();
 		end
 		function dt = get.dt(obj)
-			if isa(obj.timeDomain, "quantity.EquidistantDomain")
-				dt = obj.timeDomain.stepSize;
+			if isa(obj.domain, "quantity.EquidistantDomain")
+				dt = obj.domain.stepSize;
 			else
-				delta_t = diff(obj.timeDomain.grid);
+				delta_t = diff(obj.domain.grid);
 				assert( all( diff( delta_t ) < 1e-15 ), "Grid is not equidistant spaced" );
 				dt = delta_t(1);
 			end
@@ -39,24 +39,41 @@ classdef BasicVariable < handle
 	end
 	
 	methods ( Access = public)
-		function D = diff(obj, k)
-			% DIFF get the k-th derivative from saved derivatives
+		function D = diff(obj, domain, n)
+			% DIFF get the k-th derivative 
+			%	D = diff(obj, domain, n) tries to find the n-th derivative of obj.fun in the saved
+			%	derivatives. If it does not exist, it tries to compute it with the evaluateFunction.
+			%	If n is a non integer number, the Riemann-Liouville fractional derivative is
+			%	computed. For this, a formula is used which requires the next integer + 1 derivative
+			%	of the function to avoid the singularity under the integral. Actually, the domain is
+			%	clear, but it is used hear to be consistent with the quantity.Discrete/diff.
 			arguments
 				obj
-				k (1,1) double {mustBeInteger, mustBeNonnegative} = 0;
+				domain (1,1) quantity.Domain;
+				n (1,1) double = 0;
 			end
-						
+			
 			for l = 1:numel(obj)
-				if k == 0
+				
+				assert( domain.isequal( obj(l).domain  ),...
+					"The derivative wrt. this domain is not possible.")
+				
+				if n == 0
 					D(l) = obj(l).fun;
-				elseif k > length(obj(l).derivatives)
+					
+				elseif ~numeric.near(round(n), n)
+					% fractional derivative:
+					D(l) = obj(l).fDiff(n);
+					
+				elseif n > length(obj(l).derivatives)
 					% For some functions their maybe a possibility to compute new derivatives on the
 					% fly. For this, the evaluateFunction(k) can be used, where k is the order of
 					% the derivative.
-					D(l) = obj(l).evaluateFunction(k);
+					D(l) = obj(l).evaluateFunction(obj.domain, n);
 					obj(l).derivatives{end+1} = D(l);
 				else
-					D(l) = obj(l).derivatives{k};
+					% use the pre-computed derivatives
+					D(l) = obj(l).derivatives{n};
 				end
 			end
 			D = reshape(D, size(obj));
@@ -70,63 +87,14 @@ classdef BasicVariable < handle
 				obj,
 				k (:,1) double {mustBeInteger, mustBeNonnegative} = 0:numel(obj(1).derivatives);
 			end
-			D_ = arrayfun( @(i) obj.diff(i), k, 'UniformOutput', false);
+			D_ = arrayfun( @(i) obj.diff(obj.domain, i), k, 'UniformOutput', false);
 			
 			for i = 1:numel(k)
 				D(i,:) = D_{i};
 			end
 		end
 		
-		function D = fDiff(obj, p)
-			% FDIFF compute the p-fractional derivative
-			%	D = fDiff(obj, p) computes the fractional derivative of order p. To be
-			%	concrete the Riemann-Liouville fractional derivative with initial point t0 = lower
-			%	end of the considered time domain is computed.
-			%	Because the basisc variable must be C^infinity function, we use the form
-			%	of the fractional derivative described in eq. (2.80) in the book [Podlubny: Fractional
-			%	differential equations]. Then, integration by parts is applied to eliminate the
-			%	singularity:
-			%	
-			
-			n = ceil(p);
-			tDomain = obj(1).fun(1).domain;
-			t0 = tDomain.lower;
-			t = Discrete( tDomain );
-			tauDomain = tDomain.rename( tDomain.name + "_Hat" );
-			tau = Discrete( tauDomain );
-			if p > 0
-				
-				D = quantity.Discrete.zeros(size(obj), tDomain);
-				for k = 0:n
-					f_dk = obj.diff(k);
-					
-					if f_dk.at(t0) ~= 0
-						D = D + f_dk.at(t0) * (t - t0).^(k-p) / gamma( k-p + 1);
-					end
-					
-				end
-				assert( ~ any( isinf( obj.diff(n+1).on() ) ), "conI:BasicVariable:fDiff", ...
-					"The function has a singularity! Hence, the fractional derivative can not be computed")
-				D = D + cumInt( ( t - tau ).^(n-p) * subs(obj.diff(n+1), tDomain.name, tauDomain.name), ...
-				tauDomain, tauDomain.lower, tDomain.name) / gamma( n+1-p );
-				
-			elseif p < 0 && p > -1
-				alpha = -p;
-				D = (t - t0).^(alpha) / gamma( alpha + 1) * obj.diff(0).at(t0) + ...
-					cumInt( (t-tau).^alpha * subs(obj.diff(1), tDomain.name, tauDomain.name), ...
-					tauDomain, tauDomain.lower, tDomain.name) / gamma( alpha + 1 );
 
-			elseif p <= -1
-				alpha = -p;
-				D = cumInt( (t-tau).^(alpha-1) * subs( obj.diff(0), tDomain.name, tauDomain.name), ...
-					tauDomain, tauDomain.lower, tDomain.name) / gamma( alpha );
-				
-			elseif p == 0
-				D = [obj.fun()];	
-			else
-				error("This should not happen")
-			end
-		end
 		
 		function h = productRule(obj, phi, n)
 			% LEIBNIZ compute the k-derivative of the product of the basic variable and phi
@@ -141,7 +109,7 @@ classdef BasicVariable < handle
 			if numeric.near(round(n), n)
 				% for integer order derivatives
 				for k = 0:n
-					h = h + misc.binomial(n, k) * obj.diff(n-k) * phi.diff(t, k);
+					h = h + misc.binomial(n, k) * obj.diff(t, n-k) * phi.diff(t, k);
 				end
 			else
 				% for fractional order derivatives: [see Podlubny]
@@ -162,7 +130,7 @@ classdef BasicVariable < handle
 		
 		
 		function d = domain(obj)
-			d = obj.fun(1).domain;
+			d = obj(1).fun(1).domain;
 		end
 		
 		function plot(obj, varargin)
@@ -210,9 +178,62 @@ classdef BasicVariable < handle
 	end % methods (Access = public)
 	
 	methods (Access = protected)
-		function v = evaluateFunction(obj, k)
+		function v = evaluateFunction(obj, domain, k)
 			error('conI:signals:BasicVariable:derivative', ['Requested derivative ' num2str(k) ' is not available']);
 		end
+		
+		function D = fDiff(obj, p)
+			% FDIFF compute the p-fractional derivative
+			%	D = fDiff(obj, p) computes the fractional derivative of order p. To be
+			%	concrete the Riemann-Liouville fractional derivative with initial point t0 = lower
+			%	end of the considered time domain is computed.
+			%	Because the basisc variable must be C^infinity function, we use the form
+			%	of the fractional derivative described in eq. (2.80) in the book [Podlubny: Fractional
+			%	differential equations]. Then, integration by parts is applied to eliminate the
+			%	singularity:
+			%	
+			
+			n = ceil(p);
+			tDomain = obj(1).domain;
+			t0 = tDomain.lower;
+			t = Discrete( tDomain );
+			tauDomain = tDomain.rename( tDomain.name + "_Hat" );
+			tau = Discrete( tauDomain );
+			if p > 0
+				
+				% TODO Try to implement the polynomial part as function or as symbolic
+				D = quantity.Discrete.zeros(size(obj), tDomain);
+				for k = 0:n
+					f_dk = obj.diff(tDomain, k);
+					
+					if f_dk.at(t0) ~= 0
+						D = D + f_dk.at(t0) * (t - t0).^(k-p) / gamma( k-p + 1);
+					end
+					
+				end
+				assert( ~ any( isinf( obj.diff(tDomain, n+1).on() ) ), "conI:BasicVariable:fDiff", ...
+					"The function has a singularity! Hence, the fractional derivative can not be computed")
+				D = D + cumInt( ( t - tau ).^(n-p) * subs(obj.diff(tDomain, n+1), tDomain.name, tauDomain.name), ...
+				tauDomain, tauDomain.lower, tDomain.name) / gamma( n+1-p );
+				
+			elseif p < 0 && p > -1
+				alpha = -p;
+				D = (t - t0).^(alpha) / gamma( alpha + 1) * obj.diff(tDomain, 0).at(t0) + ...
+					cumInt( (t-tau).^alpha * subs(obj.diff(tDomain, 1), tDomain.name, tauDomain.name), ...
+					tauDomain, tauDomain.lower, tDomain.name) / gamma( alpha + 1 );
+
+			elseif p <= -1
+				alpha = -p;
+				D = cumInt( (t-tau).^(alpha-1) * subs( obj.diff(tDomain, 0), tDomain.name, tauDomain.name), ...
+					tauDomain, tauDomain.lower, tDomain.name) / gamma( alpha );
+				
+			elseif p == 0
+				D = [obj.fun()];	
+			else
+				error("This should not happen")
+			end
+		end		
+		
 	end	
 	
 end
\ No newline at end of file

+signals/GevreyFunction.m

@@ -4,7 +4,7 @@ classdef GevreyFunction < signals.BasicVariable
 	%
 	%          / 0             : t <= 0
 	%          |
-	%	f(t) = | f0 + K * h(t) : 0 < t < T
+	%   f(t) = | f0 + K * h(t) : 0 < t < T
 	%          | 
 	%          \ 0             : t >= T
 	%
@@ -19,8 +19,6 @@ classdef GevreyFunction < signals.BasicVariable
 		offset (1,1) double = 0;
 		% the number of derivatives to be considered
 		numDiff (1,1) double = 5;
-		% time domain
-		timeDomain (1,1) quantity.Domain = quantity.Domain.defaultDomain(100, "t");
 		% underlying gevery function
 		g;
 		% derivative shift
@@ -34,6 +32,10 @@ classdef GevreyFunction < signals.BasicVariable
 	
 	methods
 		function obj = GevreyFunction(optArgs)
+			% GEVREYFUNCTION initialization of a gevrey function as signals.BasicVariable
+			%	obj = GevreyFunction(optArgs) creates a gevrey function object. The properties can
+			%	be set by name-value-pairs and are "order", "gain", "offset", "numDiff",
+			%	"timeDomain", "g", "diffShift".
 			arguments
 				optArgs.order (1,1) double = 1.99;
 				optArgs.gain double = 1;
@@ -47,13 +49,14 @@ classdef GevreyFunction < signals.BasicVariable
 			% generate the underlying gevrey functions:
 			f_derivatives = cell(optArgs.numDiff + 1, 1);
 			for k = 0:optArgs.numDiff			
-				f_derivatives{k+1} = signals.GevreyFunction.evaluateFunction_p(optArgs, k);
+				f_derivatives{k+1} = ...
+					signals.GevreyFunction.evaluateFunction_p(optArgs, optArgs.timeDomain, k);
 			end
 			
 			obj@signals.BasicVariable(f_derivatives{1}, f_derivatives(2:end));	
 
 			% set the parameters to the object properties:
-			for k = fieldnames(optArgs)'
+			for k = fieldnames(rmfield(optArgs, "timeDomain"))'
 				obj.(k{:}) = optArgs.(k{:});
 			end
 			
@@ -68,7 +71,7 @@ classdef GevreyFunction < signals.BasicVariable
 	end
 	
 	methods (Static, Access = private)
-		function v = evaluateFunction_p(optArgs, k)
+		function v = evaluateFunction_p(optArgs, domain, k)
 			
 			derivativeOrder = k + optArgs.diffShift;
 			% The offset is only added to the zero-order derivative. 
@@ -76,13 +79,13 @@ classdef GevreyFunction < signals.BasicVariable
 			
 			v = quantity.Function(@(t) offset_k + optArgs.gain * ...
 					optArgs.g(t, derivativeOrder, ...
-				optArgs.timeDomain.upper, 1 ./ ( optArgs.order - 1)), optArgs.timeDomain);
+				domain.upper, 1 ./ ( optArgs.order - 1)), domain);
 		end		
 	end
 	
 	methods (Access = protected)
-		function v = evaluateFunction(obj, k)
-			v = signals.GevreyFunction.evaluateFunction_p(obj, k);
+		function v = evaluateFunction(obj, domain, k)
+			v = signals.GevreyFunction.evaluateFunction_p(obj, domain, k);
 		end
 	end
 end

+unittests/+signals/signalsTest.m

@@ -36,7 +36,7 @@ function testGevreyBump(testCase)
 % signals.
 
 g = signals.GevreyFunction('diffShift', 1);
-b = signals.gevrey.bump.dgdt_1(g.T, g.sigma, g.timeDomain.grid) * 1 / signals.gevrey.bump.dgdt_0(g.T, g.sigma, g.T);
+b = signals.gevrey.bump.dgdt_1(g.T, g.sigma, g.domain.grid) * 1 / signals.gevrey.bump.dgdt_0(g.T, g.sigma, g.T);
 
 %% Test scaling of Gevrey function
 verifyEqual(testCase, b, g.fun.on(), 'AbsTol', 1e-13);
@@ -45,7 +45,7 @@ verifyEqual(testCase, b, g.fun.on(), 'AbsTol', 1e-13);
 g = signals.GevreyFunction('diffShift', 0);
 
 orders = 1:5;
-derivatives = arrayfun(@(i) g.diff(i), orders);
+derivatives = arrayfun(@(i) g.diff(g.domain, i), orders);
 verifyEqual(testCase, [derivatives.at(0), derivatives.at(g.T)], zeros( 1, 2*length(orders)), 'AbsTol', 1e-13)
 
 end
@@ -82,29 +82,29 @@ function basicVariable_Diff_test(tc)
 	b(1,:) = signals.BasicVariable(F, {d1F, d2F});
 	b(2,:) = signals.BasicVariable(F, {d1F, d2F});
 		
-	tc.verifyEqual(b.diff(0), [F; F]);
-	tc.verifyEqual(b.diff(1).on(), on(d1F * ones(2,1)));
-	tc.verifyEqual(b.diff(2).on(), on(d2F * ones(2,1)));
-	tc.verifyError(@() b.diff(3), 'conI:signals:BasicVariable:derivative')
+	tc.verifyEqual(b.diff(T, 0), [F; F]);
+	tc.verifyEqual(b.diff(T, 1).on(), on(d1F * ones(2,1)));
+	tc.verifyEqual(b.diff(T, 2).on(), on(d2F * ones(2,1)));
+	tc.verifyError(@() b.diff(T, 3), 'conI:signals:BasicVariable:derivative')
 end
 
 function testGevreyFunction(testCase)
 	
 N = 3;
-
+tau = quantity.Domain.defaultDomain(42, "tau");
 g(1,1) = signals.GevreyFunction('order', 1.5, 'gain', -5, 'offset', 2.5, 'numDiff', N, ...
-	'timeDomain', quantity.Domain.defaultDomain(42, "tau"), 'diffShift', 0);
+	'timeDomain', tau, 'diffShift', 0);
 
 g(2,1) = signals.GevreyFunction('order', 1.9, 'gain', 5, 'offset',- 2.5, 'numDiff', N, ...
-	'timeDomain', quantity.Domain.defaultDomain(42, "tau"), 'diffShift', 0);
+	'timeDomain', tau, 'diffShift', 0);
 
-testCase.verifyEqual(g.diff(0).at(0), [2.5; -2.5])
-testCase.verifyEqual(g.diff(0).at(1), [-2.5; 2.5])
+testCase.verifyEqual(g.diff(tau, 0).at(0), [2.5; -2.5])
+testCase.verifyEqual(g.diff(tau, 0).at(1), [-2.5; 2.5])
 
 % test all the precomputed derivatives plus one higher derivative
 for k = 1:N+1
-	testCase.verifyEqual(g.diff(k).at(0), zeros(2,1))
-	testCase.verifyEqual(g.diff(k).at(1), zeros(2,1))
+	testCase.verifyEqual(g.diff(tau, k).at(0), zeros(2,1))
+	testCase.verifyEqual(g.diff(tau, k).at(1), zeros(2,1))
 end
 
 end
@@ -112,22 +112,22 @@ end
 function testGevreyFunctionTimes(testCase)
 
 N = 3;
-
+tau = quantity.Domain.defaultDomain(42, "tau");
 g(1,1) = signals.GevreyFunction('order', 1.5, 'gain', -5, 'offset', 2.5, 'numDiff', N, ...
-	'timeDomain', quantity.Domain.defaultDomain(42, "tau"), 'diffShift', 0);
+	'timeDomain', tau, 'diffShift', 0);
 
 g(2,1) = signals.GevreyFunction('order', 1.9, 'gain', 5, 'offset',- 2.5, 'numDiff', N, ...
-	'timeDomain', quantity.Domain.defaultDomain(42, "tau"), 'diffShift', 0);
+	'timeDomain', tau, 'diffShift', 0);
 
 A = [1 2; 3 4];
 
 b = A*g;
-testCase.verifyEqual(b.diff(0).at(0), A * [2.5; -2.5])
+testCase.verifyEqual(b.diff(tau, 0).at(0), A * [2.5; -2.5])
 testCase.verifyEqual(b(1).highestDerivative, g(1).highestDerivative)
 
 end
 
-function test_BasicVariable_fDiff(testCase)
+function testBasicVariableFractionalDerivative(testCase)
 	
 	t = sym('t');
 	t0 = 0;
@@ -138,7 +138,7 @@ function test_BasicVariable_fDiff(testCase)
 	F = quantity.Symbolic(( t - t0 )^beta, timeDomain);
 	b(1,:) = signals.BasicVariable(F, {F.diff("t", 1), F.diff("t", 2), F.diff("t", 3)});
 	p = 0.3;
-	testCase.verifyError( @() b.fDiff(p), "conI:BasicVariable:fDiff")
+	testCase.verifyError( @() b.diff(timeDomain, p), "conI:BasicVariable:fDiff")
 	
 	
 	%% verify regular case
@@ -150,23 +150,23 @@ function test_BasicVariable_fDiff(testCase)
 		for k = 1:5
 			p = k*alpha;
 			% compute the fractional derivative with the formula for the polynomial (see Podlubny: 2.3.4)
-			f_fdiff(k,1) = quantity.Symbolic( fractional.monomialFractionalDerivative(t, t0, beta, p), timeDomain);
+			f_diff(k,1) = quantity.Symbolic( fractional.monomialFractionalDerivative(t, t0, beta, p), timeDomain);
 			
 			% compute the fractional derivative numerically
-			b_fdiff(k,1) = b.fDiff(p);
+			b_diff(k,1) = b.diff(timeDomain, p);
 		end
-		% plot(f_fdiff - b_fdiff)
+		% plot(f_diff - b_diff)
 		
-		testCase.verifyEqual( b_fdiff.on(), f_fdiff.on(), 'AbsTol', 2e-2)
+		testCase.verifyEqual( b_diff.on(), f_diff.on(), 'AbsTol', 2e-2)
 	end
 	
 	%% verify fractional integration
-	testCase.verifyEqual( b.fDiff(0).on(), F.on() )
+	testCase.verifyEqual( b.diff(timeDomain, 0).on(), F.on() )
 	
 	for k = 1:5
 		alpha = 0.3 * k;
 		J = int( (t- sym("tau")).^(alpha -1) * subs(F.sym, "t", "tau"), "tau", t0, t) / gamma(alpha);
-		J_num = b.fDiff( -alpha );
+		J_num = b.diff( timeDomain, -alpha );
 		
 		J = quantity.Symbolic(J, timeDomain);
 		
@@ -175,7 +175,7 @@ function test_BasicVariable_fDiff(testCase)
 	
 end
 
-function test_BasicVariable_fractional_ProductRule(testCase)
+function testBasicVariableFractionalProductRule(testCase)
 
 t = quantity.Domain("t", linspace(0,1));
 % choose a polynomial basis
@@ -194,12 +194,12 @@ b = signals.BasicVariable( h(1), num2cell(h(2:end)) );
 
 for k = 1:3
 	
-	f = sym(b.diff(0)) * sym(poly(k));
+	f = sym(b.diff(t, 0)) * sym(poly(k));
 	
 	alpha = 0.3 * k;
 	n = ceil(alpha);
 	
-	frac_f = quantity.Symbolic(fractional.fDiff(f, alpha, t.lower, t.name), t);
+	frac_f = quantity.Symbolic(fractional.derivative(f, alpha, t.lower, t.name), t);
 	frac_p = b.productRule(poly(k), alpha);
 	
 	testCase.verifyEqual(frac_f.on(), frac_p.on(), 'AbsTol', 4e-3) 
@@ -215,7 +215,6 @@ phi = quantity.Symbolic( sin( sym("t") * pi), t);
 bfun = quantity.Symbolic( cos(sym("t") * pi), t);
 b = signals.BasicVariable( bfun, { bfun.diff(t, 1), bfun.diff(t, 2), bfun.diff(t, 3) });
 
-
 for k = 0:3
 	h1 = b.productRule(phi, k);
 	h2 = diff( phi * bfun, t, k);

+unittests/testFractional.m

@@ -14,7 +14,7 @@ t0 = 0.5;
 		[y, m] = fractional.monomialFractionalDerivative( t, t0, 1, p);
 		
 		% compute the fractional derivative with the function
-		f = fractional.fDiff(m, p, t0, t);
+		f = fractional.derivative(m, p, t0, t);
 		
 		% evaluate both and compare them
 		dt = 0.1;
@@ -29,7 +29,7 @@ t0 = 0.5;
 		
 		testCase.verifyEqual( Y.on(), F.on(), 'AbsTol', 1e-15);
 		
- 		h = fractional.fDiff(m, p, t0, t, 'type', 'C');
+ 		h = fractional.derivative(m, p, t0, t, 'type', 'C');
 		H = quantity.Symbolic( h, tDomain);
  		testCase.verifyEqual( H.on(), Y.on(), 'AbsTol', 1e-15);
 	end