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/+signals/BasicVariable.m b/+signals/BasicVariable.m index e3f934843c895ded16766e5a48484df921f4a1b6..2cdf83b99a2a0374af487387e5f5d0128a221044 100644 --- a/+signals/BasicVariable.m +++ b/+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 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/+signals/signalsTest.m b/+unittests/+signals/signalsTest.m index 568cee40ae6841dd6d7129a9dea04ce55e36516b..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,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); 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