diff --git a/+fractional/derivative.m b/+fractional/derivative.m
new file mode 100644
index 0000000000000000000000000000000000000000..d99fc000ac694d55dd10521b07af6373ccfe94cb
--- /dev/null
+++ b/+fractional/derivative.m
@@ -0,0 +1,39 @@
+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};
+ t0 (1,1) double;
+ t (1,1) sym;
+ optArg.type = "RL"
+end
+
+if p < 0
+ error('a must be positive! Use fintegral instead?');
+end
+
+if misc.nearInt(p)
+ h = diff(y, round(p));
+ return
+end
+
+n = ceil(p);
+tau = sym( string(t) + "_");
+
+if optArg.type == "C"
+ 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"
+ 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
+
+
+end
\ No newline at end of file
diff --git a/+fractional/fDiff.m b/+fractional/fDiff.m
deleted file mode 100644
index 4e8bceb112e82a023f634fb7c38b00e0e3c227bc..0000000000000000000000000000000000000000
--- a/+fractional/fDiff.m
+++ /dev/null
@@ -1,38 +0,0 @@
-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.
-arguments
- y sym
- p (1,1) double {mustBePositive};
- t0 (1,1) double;
- t (1,1) sym;
- optArg.type = "RL"
-end
-
-if p < 0
- error('a must be positive! Use fintegral instead?');
-end
-
-if misc.nearInt(p)
- DaY = diff(y, round(p));
- return
-end
-
-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);
- %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);
-else
- error("The derivative of the type " + optArg.type + " is not implemented. Only C for Caputo and RL for Riemann-Liouville are available")
-end
-
-
-end
\ No newline at end of file
diff --git a/+fractional/monomialFractionalDerivative.m b/+fractional/monomialFractionalDerivative.m
index bf5f5ebe61adf22c3606ce5036523cff728ce5a4..23689944550178b97f5a17d43a20dfb103a84361 100644
--- a/+fractional/monomialFractionalDerivative.m
+++ b/+fractional/monomialFractionalDerivative.m
@@ -1,7 +1,10 @@
function [y, m] = monomialFractionalDerivative(t, t0, nu, p, optArg)
%MONOMIALFRACTIONALDERIVATIVE compute fractional derivative of a monomial
% [y, m] = monomialFractionalDerivative(t, a, nu, p, optArg) compute the p-fractional derivative of
-% the monomial m = (t-a)^nu.
+% the monomial m = (t-a)^nu. At this, t must be a symbolic variable, t0 is the lower terminal as
+% double value, nu is the exponent as double, p is the fractional order as double the type of the
+% derivative can be specified by the name-value-pair "type". So far online the "RL" :
+% Riemann-Liouville derivative is implemented.
arguments
t (1,1) sym; % symbolic variable for the monomial
t0 (1,1) double; % expansion point of the monome
diff --git a/+misc/binomial.m b/+misc/binomial.m
index ce628043616125055d414ac8dec6d8c9071c1e56..a382b634772300767bf3ce5d5cb90a487d44d8ae 100644
--- a/+misc/binomial.m
+++ b/+misc/binomial.m
@@ -6,10 +6,10 @@ arguments
k (1,1) double {mustBeInteger};
end
-if p >= 0
- b = gamma(p+1) / gamma(k+1) / gamma(p-k + 1);
-else
+if p < 0 && k ~= 0
b = (-1)^k * misc.binomial( -p + k -1, k);
+else
+ b = gamma(p+1) / gamma(k+1) / gamma(p-k + 1);
end
end
diff --git a/+signals/BasicVariable.m b/+signals/BasicVariable.m
index ee5ee117729e122c7663cc25055c60d1259caf68..f735ac8a82c7081f2e561cdabd25c0f6f10d2795 100644
--- a/+signals/BasicVariable.m
+++ b/+signals/BasicVariable.m
@@ -1,5 +1,9 @@
classdef BasicVariable < handle
- %BasicVariable Class to handle basic variables for the trajectory plannning
+ %BasicVariable Class to handle quantities with a priori known derivatives
+ % This class is used to simplify the usage of quantities for which a lot of computations must be
+ % done on its derivatives. So a BasicVariable consists of a basic function and its derivatives,
+ % which are precomputed and safed. So computations can be performed without recalculation of the
+ % derivatives.
properties ( SetAccess = protected )
derivatives cell;
@@ -21,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
@@ -35,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));
@@ -66,80 +87,50 @@ 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:
- %
+
+
+ function h = productRule(obj, phi, n)
+ % LEIBNIZ compute the k-derivative of the product of the basic variable and phi
+ % h = productRule(obj, phi, n) computes the k-th derivative of obj.fun*phi, i.e.,
+ % h = d_t^n ( obj.fun(t) * phi )
+ % using the product rule. This is usefull if the derivatives of the basic variable obj
+ % are known exactly and the derivatives of phi, but not the multiplication
+ % obj.diff(0)*phi.
+ t = phi.domain;
+ h = quantity.Discrete.zeros( size(obj), t );
- 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);
+ if numeric.near(round(n), n)
+ % for integer order derivatives
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
-
+ h = h + misc.binomial(n, k) * obj.diff(t, n-k) * phi.diff(t, k);
end
-
- 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")
+ % for fractional order derivatives: [see Podlubny]
+ % p \~~ oo / p \ (k) p-k
+ % D (phi(t) f(t)) = > | | | phi (t) D f(t)
+ % a t /__ k=0 \ k / a t
+ h = misc.binomial( n, 0 ) * phi * obj.fDiff(n);
+ h0 = h;
+ k = 1;
+ % do the iteration on the loop until the difference is very small.
+ while max( abs( h0 ) ) / max( abs( h ) ) > 1e-9
+ h0 = misc.binomial(n,k) * phi.diff(t, k) * obj.fDiff( n - k);
+ h = h + h0;
+ k = k+1;
+ end
end
end
- function D = fLeibniz(obj, phi, p)
- % FLEIBNIZ compute the fractional derivative of multiplication of two functions:
- % D^p ( phi(t) * obj.fun(t) )
- % = sum_i=0^infty (p, k) * dt^k phi(t) * D^(p-k) obj.fun(t)
- %
-
- tau = phi.domain.name;
- N = 10;
- D = misc.binomial( p, 0 ) * phi * obj.fDiff(p);
- for k = 1:N % TODO@ferdinand: make while loop with lower threshold.
- D = D + misc.binomial(p,k) * phi.diff(tau, k) * obj.fDiff( p - k);
- end
- end
-
function d = domain(obj)
- d = obj.fun(1).domain;
+ d = obj(1).fun(1).domain;
end
function plot(obj, varargin)
@@ -187,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
+end
diff --git a/+signals/GevreyFunction.m b/+signals/GevreyFunction.m
index c817447e19dbd5b0dca9c48d357f1c31b12b11be..49ae6323fd7f7d41f40fcd8537a727b79019d96e 100644
--- a/+signals/GevreyFunction.m
+++ b/+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
diff --git a/+unittests/+misc/testMisc.m b/+unittests/+misc/testMisc.m
index 5a18830a101d09635364304f238988dd7a7416f2..afd1974b849dce5f59031fd0167387e6b41d66dc 100644
--- a/+unittests/+misc/testMisc.m
+++ b/+unittests/+misc/testMisc.m
@@ -19,7 +19,7 @@ end
testCase.verifyEqual( misc.binomial(2.5, 2), 1.875)
testCase.verifyEqual( misc.binomial(-1, 4), (-1)^4)
-for k = 1:3
+for k = -1:3
% test non integer over zero
testCase.verifyEqual( misc.binomial( k*0.3, 0), 1);
% test non integer over one
diff --git a/+unittests/+signals/signalsTest.m b/+unittests/+signals/signalsTest.m
index 26d892466421e96fd3d76960828e710c1fe1e0d2..6980acce85b50c7b60eafcf1d046b302ef2effab 100644
--- a/+unittests/+signals/signalsTest.m
+++ b/+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,89 +112,136 @@ 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;
timeDomain = quantity.Domain('t', linspace(t0, t0+1, 1e2));
- beta = 2;
-
- F = quantity.Symbolic(( t - t0 )^beta, timeDomain);
- d1F = F.diff("t", 1);
- d2F = F.diff("t", 2);
+ %% verify singularity case
+ beta = 1.5;
+ 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;
- nu = 2;
- for k = 1:5
- % compute the fractional derivative with the formula for the polynomial (see Podlubny: 2.3.4)
- f_fdiff(k,1) = quantity.Symbolic( gamma(1+nu) / gamma( 1 + nu - k*p ) * ( t - t0 ).^(nu - k*p), timeDomain);
-
- % compute the fractional derivative numerically
- b_fdiff(k,1) = b.fDiff(k*p);
- end
- % plot(f_fdiff - b_fdiff)
+ testCase.verifyError( @() b.diff(timeDomain, p), "conI:BasicVariable:fDiff")
- testCase.verifyEqual( b_fdiff.on(), f_fdiff.on(), 'AbsTol', 1e-2)
+ %% verify regular case
+ for beta = [1, 2, 3]
+ F = quantity.Symbolic(( t - t0 )^beta, timeDomain);
+ b(1,:) = signals.BasicVariable(F, {F.diff("t", 1), F.diff("t", 2), F.diff("t", 3)});
+
+ alpha = 0.3;
+ for k = 1:5
+ p = k*alpha;
+ % compute the fractional derivative with the formula for the polynomial (see Podlubny: 2.3.4)
+ f_diff(k,1) = quantity.Symbolic( fractional.monomialFractionalDerivative(t, t0, beta, p), timeDomain);
+
+ % compute the fractional derivative numerically
+ b_diff(k,1) = b.diff(timeDomain, p);
+ end
+ % plot(f_diff - b_diff)
+
+ 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);
- testCase.verifyEqual( J.on() , J_num.on, 'AbsTol', 2e-3 );
+ testCase.verifyEqual( J.on() , J_num.on, 'AbsTol', 3e-3 );
end
end
-function test_BasicVariable_fLeibniz(testCase)
+function testBasicVariableFractionalProductRule(testCase)
t = quantity.Domain("t", linspace(0,1));
% choose a polynomial basis
N = 5;
-for k = 1:N
- poly(k) = quantity.Symbolic(sym("t").^(k-1), t);
+
+
+beta = 1;
+h(1,1) = quantity.Symbolic(sym("t").^beta, t);
+poly(1,1) = quantity.Symbolic( sym("t"), t);
+for k = 2:N
+ h(k,1) = h(k-1).diff();
+ poly(k,1) = quantity.Symbolic( sym("t")^(k-1) / factorial(k-1), t);
end
-p = signals.BasicVariable( poly(2), num2cell(poly(3:end)) );
+b = signals.BasicVariable( h(1), num2cell(h(2:end)) );
for k = 1:3
- f = sym(p.diff(0)) * sym(poly(1));
+ f = sym(b.diff(t, 0)) * sym(poly(k));
alpha = 0.3 * k;
n = ceil(alpha);
+
+ 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)
+end
- frac_f = diff( ...
- int( (sym("t") - sym("tau"))^(n-alpha-1) * subs(f, "t", "tau"), "tau", 0, "t"), ...
- "t", n) / gamma(n - alpha);
+end
- frac_f = quantity.Symbolic(frac_f, t);
-
- frac_p = p.fLeibniz(poly(1), alpha);
+
+function test_BasicVariable_ProductRule(testCase)
+
+t = quantity.Domain("t", linspace(0,1,11));
+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);
+ testCase.verifyEqual( h1.on, h2.on, 'AbsTol', 1e-13 );
end
end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/+unittests/testFractional.m b/+unittests/testFractional.m
index db6aff94c7886670c45c02976cb6f5b873b9cb26..1ebc3e1da83040c3fa8d2aeeed35406fa52c0474 100644
--- a/+unittests/testFractional.m
+++ b/+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