Commit e0d4c215 authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

simplification of quantity.Symbolic/int and /cumInt

it is now adopted to the implementation of the quantity.Discrete/int/cumInt
parent 3369855c
function F = cumIntegral(f, grid)
% CUMINTEGRAL cumulative integral for function handles
% F = cumIntegral(f, lower, upper)
arguments
f function_handle;
grid (:,1) double;
end
F = zeros(numel(grid), 1);
for i = 2:numel(grid)
F(i) = F(i-1) + integral(f, grid(i-1), grid(i));
end
end
......@@ -16,18 +16,6 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
% number of discretization points == gridLength
n (1,1) double;
% Definition of the end-points of the considered domain. By the
% logical array endpointsIncluded it can be specified if the upper
% and lower endpoint of the domain should be included or not. To be
% specific
% D = (a, b) or D = [a, b]
% or are mixture, can be specified by the logical values. TRUE
% represents inclusion, i.e., "[" or "]" and FALSE represents
% exclusion, i.e., "(" or ")".
% EXAMPLE: to specify the bounds D = (a, b], the property must be
% set as
% endpointsIncluded = [false, true]
endpointsIncluded (2,1) logical;
end
properties (Dependent)
......@@ -40,20 +28,18 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
end
methods
function obj = Domain(name, grid, optArg)
function obj = Domain(name, grid)
%DOMAIN initialize the domain
arguments
name string = "";
grid double = [];
optArg.endpointsIncluded (2,1) logical = [true true];
end
if nargin >= 1
obj.grid = grid;
obj.name = name;
obj.n = length(obj.grid);
obj.endpointsIncluded = optArg.endpointsIncluded;
end
end % Domain()
......@@ -75,6 +61,9 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
d = quantity.Domain(obj.name, obj.grid);
end
function s = string(obj)
s = obj.name;
end
end
methods (Access = public)
......
......@@ -71,6 +71,16 @@ classdef Function < quantity.Discrete
end
end
function Q = quantity.Discrete(obj, varargin)
myParser = misc.Parser();
myParser.addParameter('domain', obj(1).domain);
myParser.addParameter('name', obj(1).name);
myParser.parse(varargin{:});
Q = quantity.Discrete(obj.on(), ...
'domain', myParser.Results.domain, ...
'name', myParser.Results.name);
end
%-----------------------------
% --- overloaded functions ---
%-----------------------------
......@@ -135,11 +145,41 @@ classdef Function < quantity.Discrete
end
function I = int(S, a, b)
I = zeros(size(S));
for k = 1:numel(S)
I(k) = integral(@(varargin)S(k).evaluateFunction(varargin{:}), a, b, 'AbsTol', 1e-9);
function I = cumInt(obj, domain, lowerBound, upperBound)
% input parser since some default values depend on intGridName.
myParser = misc.Parser;
myParser.addRequired('domain', @(d) obj(1).domain.gridIndex(d) ~= 0);
myParser.addRequired('lowerBound', ...
@(l) isnumeric(l) || ischar(l) || isstring(l));
myParser.addRequired('upperBound', ...
@(l) isnumeric(l) || ischar(l) || isstring(l));
myParser.parse(domain, lowerBound, upperBound)
% get the desired domain. This step is necessary, as the input
% argument DOMAIN can be either a string or a quantity.Domain
% object.
domain = obj(1).domain.find(domain);
for k = 1:numel(obj)
I_val = numeric.cumIntegral( ...
@(varargin)obj(k).evaluateFunction(varargin{:}), ...
domain.grid);
result = quantity.Discrete(I_val, 'domain', obj(1).domain);
% int_lowerBound^upperBound f(.) =
% F(upperBound) - F(lowerBound)
I(k) = result.subs(domain, upperBound) ...
- result.subs(domain, lowerBound);
if isa(result, 'quantity.Discrete')
result.setName("int(" + obj(1).name + ")");
end
end
I = reshape(I, size(obj));
end
end
......
......@@ -174,17 +174,17 @@ classdef Symbolic < quantity.Function
res = ~(A==B);
end
function Q = quantity.Discrete(obj, varargin)
% Cast of a quantity.Symbolic object into a quantity.Discrete
% object.
myParser = misc.Parser();
myParser.addParameter('domain', obj(1).domain);
myParser.addParameter('name', obj(1).name);
myParser.parse(varargin{:});
Q = quantity.Discrete(obj.on(), ...
'domain', myParser.Results.domain, ...
'name', myParser.Results.name);
end
% function Q = quantity.Discrete(obj, varargin)
% % Cast of a quantity.Symbolic object into a quantity.Discrete
% % object.
% myParser = misc.Parser();
% myParser.addParameter('domain', obj(1).domain);
% myParser.addParameter('name', obj(1).name);
% myParser.parse(varargin{:});
% Q = quantity.Discrete(obj.on(), ...
% 'domain', myParser.Results.domain, ...
% 'name', myParser.Results.name);
% end
function f = function_handle(obj)
f = matlabFunction(sym(obj));
end
......@@ -721,80 +721,24 @@ classdef Symbolic < quantity.Function
b = obj.variable;
z = num2cell( obj(1).variable );
end
if ~iscell(z)
z = {z};
end
if obj.nargin == numel(z) && isnumeric(a) && isnumeric(b)
C = int@quantity.Function(obj, a, b);
% for k = 1:numel(obj)
% C(k) = double( int( sym(obj(k)), z{:}, a, b) );
% end
% C = reshape(C, size(C));
else
% compute the symbolic integration
I = int(obj.sym, z{:}, a, b);
variableI = quantity.Symbolic.getVariable(I);
% avoid empty variableI:
if ischar(a) && ~any(strcmp(string(variableI(:)), a))
variableI = [variableI(:); sym(a)];
end
if ischar(b) && ~any(strcmp(string(variableI(:)), b))
variableI = [variableI(:); sym(b)];
end
% get gridName from variableI
if iscell(variableI)
gridNameI = cellfun(@(v) char(v), variableI);
else%if isvector(variableI)
gridNameI = arrayfun(@(v) char(v), variableI, 'UniformOutput', false);
end
% create grid for result
gridI = cell(1, numel(gridNameI));
for it = 1 : numel(gridNameI)
gridIdx = obj(1).domain.gridIndex(gridNameI{it});
if (gridIdx ~= 0) && ~isempty(gridIdx)
gridI{it} = obj.gridOf(gridNameI{it});
else
% add default grid
gridI{it} = linspace(0, 1, 101);
end
end
try
% if the integral could not be calculated explicitly, the
% integral() function for functions is called. This can only deal
% with numeric integral bounds. Therefore, this try-catch block
% is needed.
C = quantity.Symbolic(I, ...
'grid', gridI, 'gridName', gridNameI, ...
'name', "int(" + obj(1).name, ")");
catch
C = cumInt(quantity.Discrete(obj), z, a, b);
end
% obj.copy();
% for k = 1:numel(obj)
% C(k).valueSymbolic = I(k);
% C(k).valueDiscrete = [];
% C(k).valueContinuous = obj.setValueContinuous(I(k));
% end
end
C = obj.cumInt(z, a, b);
end % int()
function result = cumInt(obj, domain, lowerBound, upperBound)
function C = cumInt(obj, domain, lowerBound, upperBound)
% CUMINT cumulative integration
% result = cumInt(obj, domain, lowerBound, upperBound)
% performes the integration over 'obj' for the 'domain' and
% the specified 'bounds'. 'domain' must be a gridName of the
% object. 'lowerBound' and 'upperBound' define the boundaries
% of the integration domain. These can be either doubles or
% gridNames. If it is double, the bound is constant. If it is
% a gridName, the bound is variable and the result is a
% function dependent of this variable.
% Since cumInt can be done with the usual Symbolic int(), this method
% just checks the input and calls int().
% result = cumInt(obj, domain, lowerBound, upperBound)
% performes the integration over 'obj' for the 'domain' and
% the specified 'bounds'. 'domain' must be a gridName of the
% object. 'lowerBound' and 'upperBound' define the boundaries
% of the integration domain. These can be either doubles or
% gridNames. If it is double, the bound is constant. If it is
% a gridName, the bound is variable and the result is a
% function dependent of this variable. Since cumInt can be
% done with the usual Symbolic int(), this method just checks
% the input and calls int().
% input parser since some default values depend on intGridName.
%input parser since some default values depend on intGridName.
myParser = misc.Parser;
myParser.addRequired('domain', @(d) obj(1).domain.gridIndex(d) ~= 0);
myParser.addRequired('lowerBound', ...
......@@ -803,10 +747,42 @@ classdef Symbolic < quantity.Function
@(l) isnumeric(l) || ischar(l) || isstring(l) );
myParser.parse(domain, lowerBound, upperBound)
% call int to perform integration
result = int(obj, {domain}, lowerBound, upperBound);
% get the desired domain. This step is necessary, as the input
% argument DOMAIN can be either a string or a quantity.Domain
% object.
domain = obj(1).domain.find(domain);
% do the integration element-wise
if obj.nargin == numel(domain) && isnumeric(lowerBound) && isnumeric(upperBound)
% if the integration has fixed upper and lower bounds, and
% should be applied for all independent variables, use the
% implementation based on function handles, because it is
% much faster:
C = cumInt@quantity.Function(obj, domain, lowerBound, upperBound);
else
try
% try to compute the symbolic integration
I = quantity.Symbolic(int(obj.sym, domain.name), 'domain', obj(1).domain);
C = subs(I, domain.name, upperBound) - ...
subs(I, domain.name, lowerBound);
catch ex
try
% try to compute the integration by function
% handles
C = cumInt@quantity.Function(obj, domain, lowerBound, upperBound);
catch ex
% last fall back to compute the integration by
% discrete values. This should work always
C = cumInt(quantity.Discrete(obj), domain, lowerBound, upperBound);
end
end
end
end
end
methods (Access = protected)
......
function [tests] = nearTest()
tests = functiontests(localfunctions);
end
function testCumIntegral(testCase)
Z = linspace(0,1)';
f = @(z) sin(z * pi);
F = numeric.cumIntegral(f, Z);
F_sym = int( sym(f) );
F_eval = double( subs( F_sym , 'z', Z ) ) - double( subs( F_sym, 'z', 0) );
testCase.verifyEqual(F, F_eval, 'AbsTol', 3e-16);
end
\ No newline at end of file
......@@ -4,18 +4,39 @@ function [tests] = testFunction()
tests = functiontests(localfunctions);
end
function testCastSymbolic2Function(testCase)
function testInt(testCase)
z = quantity.Domain("z", linspace(0,1));
f = quantity.Function( @(z) sin( sin( z ) ), 'domain', z);
F_fn = int(f, z, 0, "z");
F_dk = int(quantity.Discrete(f), z, 0, "z");
testCase.verifyEqual( F_fn.on, F_dk.on, 'AbsTol', 6e-6)
f2 = quantity.Function({@(z)sin(z); @(z)cos(z)}, 'domain', z);
F2_fn = int(f2, z, "zeta", "z");
F2_dk = int( quantity.Discrete(f2), z, "zeta", "z");
testCase.verifyEqual( F2_fn.on, F2_dk.on, 'AbsTol', 7.5e-6);
end
function testCastFunction2Discrete(testCase)
z = linspace(0, 2*pi)';
Z = sym("z");
s = quantity.Symbolic([sin(Z); cos(Z)], 'grid', z, 'gridName', 'z');
f = quantity.Function(s);
f = quantity.Function({@(z)sin(z); @(z)cos(z)}, 'grid', z, 'gridName', 'z');
d = quantity.Discrete(f);
testCase.verifyTrue(all(size(f) == size(s)));
testCase.verifyTrue(all(size(d) == size(f)));
end
function testTimesSymbolic(testCase)
z = linspace(0, 2*pi)';
......
......@@ -16,6 +16,18 @@ myVectorSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1; 0.9-z/2], ...
testCase.verifyTrue( myVectorSymbolic.vec2diag.near(myMatrixSymbolic) );
end
function testCastSymbolic2Function(testCase)
z = linspace(0, 2*pi)';
Z = sym("z");
s = quantity.Symbolic([sin(Z); cos(Z)], 'grid', z, 'gridName', 'z');
f = quantity.Function(s);
testCase.verifyTrue(all(size(f) == size(s)));
end
function testSymbolicEvaluation(tc)
syms z
myGrid = linspace(0, 1, 7);
......@@ -107,6 +119,16 @@ fCumInt2Cum = cumInt(integrandSymbolic, 's', tGrid(1), 't') ...
testCase.verifyEqual(fCumInt2Bcs.on(fCumInt2Bcs(1).grid, fCumInt2Bcs(1).gridName), ...
fCumInt2Cum.on(fCumInt2Bcs(1).grid, fCumInt2Bcs(1).gridName), 'AbsTol', 100*eps);
%% test not an integrable function
f = quantity.Symbolic([ sin( sin( s ) ); cos( cos( t ) )], ...
'grid', {tGrid, sGrid}, 'gridName', {'t', 's'});
F_sym = cumInt(f, 's', 0, 1);
F_dic = cumInt( quantity.Discrete(f), 's', 0, 1);
testCase.verifyEqual( F_sym.on(), F_dic.on())
end
function testScalarPlusMinusQuantity(testCase)
......@@ -455,7 +477,7 @@ end
end
function testCastToQuantity(testCase)
function testCastSymbolic2Discrete(testCase)
syms z t
x = quantity.Symbolic(sin(z * t * pi), ...
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment