Commit b9eb20e8 authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

quantity.Discrete.diff: now supports non-equidistant grids

parent 9d26d8c5
......@@ -2415,25 +2415,27 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
function result = diff_inner(obj, k, diffGridName)
gridSelectionIndex = obj(1).domain.gridIndex(diffGridName);
gridIndex = obj(1).domain.gridIndex(diffGridName);
spacing = gradient(obj(1).grid{gridSelectionIndex}, 1);
assert(numeric.near(spacing, spacing(1)), ...
'diff is currently only implemented for equidistant grid');
permutationVector = 1 : (numel(obj(1).grid)+ndims(obj));
objDiscrete = permute(obj.on(), ...
[permutationVector(gridSelectionIndex), ...
permutationVector(permutationVector ~= gridSelectionIndex)]);
[permutationVector(gridIndex), ...
permutationVector(permutationVector ~= gridIndex)]);
[spacing, idx] = getSpacing(obj);
if size(objDiscrete, 2) == 1
derivativeDiscrete = gradient(objDiscrete, spacing(1));
if iscolumn(objDiscrete)
derivativeDiscrete = gradient(objDiscrete, ...
spacing{idx == gridIndex}, ...
spacing{idx ~= gridIndex});
else
[~, derivativeDiscrete] = gradient(objDiscrete, spacing(1));
spacing = [spacing(idx == gridIndex); spacing(idx ~= gridIndex)];
[~, derivativeDiscrete] = gradient(objDiscrete, ...
spacing{2}, spacing{idx~=2});
end
rePermutationVector = [2:(gridSelectionIndex), ...
1, (gridSelectionIndex+1):ndims(derivativeDiscrete)];
rePermutationVector = [2:(gridIndex), ...
1, (gridIndex+1):ndims(derivativeDiscrete)];
result = quantity.Discrete(...
permute(derivativeDiscrete, rePermutationVector), ...
'size', size(obj), 'domain', obj(1).domain, ...
......@@ -2444,7 +2446,27 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% recursivly until the first-order derivative is reached
result = result.diff(diffGridName, k-1);
end
end
end % diff_inner()
function [mySpace, idx] = getSpacing(obj)
% getSpacing returns a cell array of the spacing needed for gradient(),
% in diff_inner including a spacing for the array dimensions of obj.
mySpace = cell(numel(obj(1).domain), 1);
for it = 1 : numel(obj(1).domain)
mySpace{it, 1} = obj(1).domain(it).grid;
end
if isscalar(obj)
% do nothing
elseif iscolumn(obj)
mySpace{it+1, 1} = (1 : 1 : numel(obj)).';
else % ismatrix(obj)
for jt = 1 : ndims(obj)
mySpace{it+jt, 1} = (1 : 1 : size(obj, jt)).';
end
end
idx = (1 : 1 : numel(mySpace)).';
end % getSpacing(obj)
function [idx, permuteGrid] = computePermutationVectors(a, b)
% Computes the required permutation vectors to use
......
......@@ -641,10 +641,20 @@ testCase.verifyEqual(ab.on(), ones(11, 1));
testCase.verifyEqual(ac.on(), ones(21, 1));
end
function testDiffNonequidistant(tc)
z = quantity.Domain("z", linspace(0, 1, 11).^2);
quan1 = quantity.Discrete(ones(z.n, 1), "domain", z);
quanZsqrt = quantity.Discrete(1+sqrt(z.grid), "domain", z);
quanZZ = quantity.Discrete(1+z.grid, "domain", z);
tc.verifyEqual(quan1.diff.on(), 0*z.grid);
tc.verifyEqual(quanZZ.diff.on(), 1+0*z.grid, 'AbsTol', 10*eps);
end % testDiffNonequidistant()
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]);
'gridName', 'z', 'size', [2, 1]);
constantDiff = diff(constant);
testCase.verifyEqual(constantDiff.on(), [zeros(11, 1), ones(11, 1)], 'AbsTol', 10*eps);
......
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