1% Sparse operator for differentiation. 2% 3% diffx = sparse_diff(sz1,m,cyclic) 4% 5% Sparse operator for differentiation along dimension m for "collapsed" matrix 6% of the size sz1. 7% 8% Input: 9% sz1: size of rhs 10% m: dimension to differentiate 11% cyclic: true if domain is cyclic along dimension m. False is the 12% default value 13 14function diffx = sparse_diff(sz1,m,cyclic) 15 16if nargin == 2 17 cyclic = false; 18end 19 20n1 = prod(sz1); 21sz2 = sz1; 22 23if ~cyclic 24 sz2(m) = sz2(m)-1; 25end 26 27n2 = prod(sz2); 28 29n = length(sz1); 30 31vi = cell(1,n); 32for i=1:n 33 vi{i} = 1:sz2(i); 34end 35 36IJ = cell(n,1); 37 38[IJ{:}] = ndgrid(vi{:}); 39 40for i=1:n 41 IJ{i} = IJ{i}(:); 42end 43 44L1 = [1:n2]'; 45 46L2 = sub2ind(sz1,IJ{:}); 47one = ones(size(L1)); 48 49IJ{m} = IJ{m} + 1; 50 51if cyclic 52 IJ{m} = mod(IJ{m}-1,sz1(m))+1; 53end 54 55 56L2o = sub2ind(sz1,IJ{:}); 57 58diffx = sparse( ... 59 [L1; L1; ], ... 60 [L2; L2o; ], ... 61 [-one; one ], n2 , n1 ); 62 63 64 65% Copyright (C) 2009,2012 Alexander Barth <a.barth@ulg.ac.be> 66% 67% This program is free software; you can redistribute it and/or modify 68% it under the terms of the GNU General Public License as published by 69% the Free Software Foundation; either version 2 of the License, or 70% (at your option) any later version. 71% 72% This program is distributed in the hope that it will be useful, 73% but WITHOUT ANY WARRANTY; without even the implied warranty of 74% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 75% GNU General Public License for more details. 76% 77% You should have received a copy of the GNU General Public License 78% along with this program; If not, see <http://www.gnu.org/licenses/>. 79 80 81% LocalWords: diffx sz 82