Commit 606ec16c authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

Allow the grid to be changed allow differentiation of constant values

parent 01053b78
......@@ -3,14 +3,6 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% Discrete evaluation of the continuous quantity
valueDiscrete double;
% Grid for the evaluation of the continuous quantity. For the
% example with the function f(x,t), the grid would be
% {[<spatial domain>], [<temporal domain>]}
% whereas <spatial domain> is the discret description of the
% spatial domain and <temporal domain> the discrete description of
% the temporal domain.
grid; % in set.grid it is ensured that, grid is a (1,:)-cell-array
% Name of the domains that generate the grid.
gridName {mustBe.unique};
......@@ -24,6 +16,15 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
properties
% Grid for the evaluation of the continuous quantity. For the
% example with the function f(x,t), the grid would be
% {[<spatial domain>], [<temporal domain>]}
% whereas <spatial domain> is the discret description of the
% spatial domain and <temporal domain> the discrete description of
% the temporal domain.
grid; % in set.grid it is ensured that, grid is a (1,:)-cell-array
% ID of the figure handle in which the handle is plotted
figureID double = 1;
......@@ -478,10 +479,21 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% the basic values for the initialization of new
% quantity.Discrete values
isAquantityDiscrete = cellfun(@(o) isa(o, 'quantity.Discrete'), objCell);
objIdx = find(isAquantityDiscrete, 1);
obj = objCell{objIdx};
isEmpty = cellfun(@(o) isempty(o), objCell);
objIdx = find(isAquantityDiscrete & (~isEmpty), 1);
for k = 1:numel(objCell)
if all(isEmpty)
% if there are only empty entries, nothing can be
% concatenated, so a new empty object is initialized.
s = cellfun(@(o) size(o), objCell, 'UniformOutput', false);
S = sum(cat(3, s{:}), 3);
c = quantity.Discrete.empty(S);
return
else
obj = objCell{objIdx};
end
for k = 1:numel(objCell(~isEmpty))
if isa(objCell{k}, 'quantity.Discrete')
o = objCell{k};
......@@ -504,10 +516,10 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
end
[fineGrid, fineGridName] = getFinestGrid(objCell{:});
[fineGrid, fineGridName] = getFinestGrid(objCell{~isEmpty});
for it = 1 : (numel(varargin) + 1) % +1 because the first entry is a
objCopyTemp = copy(objCell{it});
objCell{it} = objCopyTemp.changeGrid(fineGrid, fineGridName);
objCell{it} = objCell{it}.changeGrid(fineGrid, fineGridName);
end
assertSameGrid(objCell{:});
argin = [{dim}, objCell(:)'];
......@@ -654,7 +666,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if nargin == 1 || isempty(gridName2Replace)
% if gridName2Replace is empty, then nothing must be done.
solution = obj;
elseif isempty(obj);
elseif isempty(obj)
% if the object is empty, nothing must be done.
solution = obj;
else
......@@ -967,6 +979,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% gridName in the obj properties remains unchanged, only the
% data points are exchanged.
if isempty(obj)
newObj = obj.copy();
return;
end
gridIndexNew = obj(1).gridIndex(gridNameNew);
......@@ -1327,18 +1340,38 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
function result = diff(obj, gridName)
% diff applies the derivative for the variable specified with
function result = diff(obj, k, gridName)
% diff applies the kth-derivative for the variable specified with
% the input gridName to the obj.
if nargin == 1
gridName = obj(1).gridName{1};
k = 1;
end
if obj.isConstant && isempty(obj(1).gridName)
result = quantity.Discrete(zeros(size(obj)), ...
'size', size(obj), 'grid', obj(1).grid, ...
'gridName', obj(1).gridName, ...
'name', ['(d_{.}', obj(1).name, ')']);
return
end
if nargin <= 2
gridName = obj(1).gridName{:};
end
% if a higher order derivative is requested, call the function
% recursivly until the first-order derivative is reached
if k > 1
assert(isnumeric(k))
obj = obj.diff(k-1, gridName);
end
gridSelector = strcmp(obj(1).gridName, gridName);
gridSelectionIndex = find(gridSelector);
spacing = gradient(obj(1).grid{gridSelectionIndex}, 1);
assert(numeric.near(spacing, spacing(1)), ...
'diff is currently only implemented for equidistand grid');
'diff is currently only implemented for equidistant grid');
permutationVector = 1 : (numel(obj(1).grid)+ndims(obj));
objDiscrete = permute(obj.on(), ...
......@@ -1478,14 +1511,15 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
[gridJoinedLength(0~=bGridIndex), gridJoinedLength(~bGridIndex), ...
size(B)]);
% % permute bDiscrete such that grids are in the order specified
% % by gridNameJoined.
[bGridNameSorted, bSortIndex] = sort([gridNameJoined(bGridIndex > 0), gridNameJoined(~bGridIndex)]);
bDiscrete = permute(bDiscrete, [bSortIndex, numel(bSortIndex)+(1:ndims(B))]);
% % permute bDiscrete such that grids are in the order specified
% % by gridNameJoined.
[bGridNameSorted, bSortIndex] = sort([gridNameJoined(bGridIndex > 0), gridNameJoined(~bGridIndex)]);
bDiscrete = permute(bDiscrete, [bSortIndex, numel(bSortIndex)+(1:ndims(B))]);
assert(isequal(aGridNameSorted, bGridNameSorted) ...
&& isequal(aGridNameSorted, gridNameJoined), ...
'Sorting grids went wrong.')
% create result object
C = quantity.Discrete(aDiscrete + bDiscrete, ...
'grid', gridJoined, 'gridName', gridNameJoined, ...
......@@ -1498,6 +1532,10 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
[C.name] = deal([A(1).name, '-', B(1).name]);
end
function C = uminus(A)
C = (-1) * A;
end
function [P, supremum] = relativeErrorSupremum(A, B)
assert(numel(A) == numel(B), 'Not implemented')
......
......@@ -223,7 +223,7 @@ quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 51)}, .
'gridName', 'z', 'name', 'bDiscrete', 'size', size(quanBSym));
solutionBDiscrete = quanBDiscrete.solveDVariableEqualQuantity();
myGrid = solutionBDiscrete.grid{1};
solutionBDiscreteDiff = solutionBDiscrete.diff('s');
solutionBDiscreteDiff = solutionBDiscrete.diff(1, 's');
quanOfSolutionOfS = zeros(length(myGrid), length(myGrid), 1);
for icIdx = 1 : length(myGrid)
quanOfSolutionOfS(:, icIdx, :) = quanBSym.on(solutionBDiscrete.on({myGrid, myGrid(icIdx)}));
......@@ -254,8 +254,22 @@ function testDiff1d(testCase)
%% 1d
constant = quantity.Discrete([2*ones(11, 1), linspace(0, 1, 11).'], 'grid', linspace(0, 1, 11), ...
'gridName', 'z', 'name', 'constant', 'size', [2, 1]);
constantDiff = diff(constant);
testCase.verifyEqual(constantDiff.on(), [zeros(11, 1), ones(11, 1)], 'AbsTol', 10*eps);
z = linspace(0,pi)';
sinfun = quantity.Discrete(sin(z), 'grid', z, 'gridName', 'z');
% do the comparison on a smaller grid, because the numerical derivative is
% very bad a the boundarys of the domain.
Z = linspace(0.1, pi-0.1)';
testCase.verifyTrue(numeric.near(sinfun.diff().on(Z), cos(Z), 1e-3));
testCase.verifyTrue(numeric.near(sinfun.diff(2).on(Z), -sin(Z), 1e-3));
testCase.verifyTrue(numeric.near(sinfun.diff(3).on(Z), -cos(Z), 1e-3));
end
function testDiffConstant2d(testCase)
......@@ -264,9 +278,9 @@ function testDiffConstant2d(testCase)
myQuantity = quantity.Discrete(cat(3, 2*ones(11, 21), zNdgrid, zetaNdgrid), ...
'grid', {linspace(0, 1, 11), linspace(0, 1, 21)}, ...
'gridName', {'z', 'zeta'}, 'name', 'constant', 'size', [3, 1]);
myQuantityDz = diff(myQuantity, 'z');
myQuantityDzeta = diff(myQuantity, 'zeta');
myQuantityDzeta2 = diff(myQuantity);
myQuantityDz = diff(myQuantity, 1, 'z');
myQuantityDzeta = diff(myQuantity, 1, 'zeta');
myQuantityDzeta2 = diff(myQuantity, 1);
% constant
testCase.verifyEqual(myQuantityDz(1).on(), zeros(11, 21));
......
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