# gradientNeumann

Discretization of the gradient using standard finite differences with Neumann boundary conditions.

## Contents

- Syntax
- Description
- Example
- Input Arguments
- Output Arguments
- More About
- References
- See Also
- Code

## Syntax

    res = gradientNeumann(xnRVD,h,N,G,seti)

## Description

`res = gradientNeumann(xnRVD,h,N,G,seti)` computes the discretized gradient `res` of `xnRVD` (matrix on down scaled ROI stored as a vector) with length of the infinitesimal element `h` and discretization points for each dimension `N`.

## Example

From setFuncsPda.html:

gradientNeumann(xnRVD,seti.hInv,seti.nInv,seti.GInv,seti)

## Input Arguments

- `xnRVD` : matrix on down scaled ROI stored as a real vector (i.e. in the code: firstly we identify with and secondly reshape the vector to a matrix in the down scaled ROI).
- h : length of the infinitesimal element, see setGrid.html
- N : discretization points for each dimension, see setGrid.html
- G : function to reshape a vector to a matrix, see setReshapeVecMat.html
- seti : structural array, in this file we need the transformation operators `seti.T` and `seti.S` to identify with .

## Output Arguments

- res : gradient of xnRVD (real array of size 2*2 x N x N in case of 2D and size 2*3 x N x N x N in case of 3D).

## More About

Discretization of using standard finite differences with Neumann boundary conditions, see [1, Sec. 6.1] or [2, Sec. 4.3].

We present the 2D case (this can be extended to 3D).

with

- .

and

- .

## References

- [1] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120-145, 2011. URL: https://doi.org/10.1007/s10851-010-0251-1.
- [2] Florian Bürgel, Kamil S. Kazimierski, and Armin Lechleiter. A sparsity regularization and total variation based computational framework for the inverse medium problem in scattering. Journal of Computational Physics, 339:1-30, 2017. URL: https://doi.org/10.1016/j.jcp.2017.03.011.

## See Also

- gradientNeumannAdj.html

## Code

    function res = gradientNeumann(xnRVD,h,N,G,seti)

    % xnRVD: real, stored as vector, down scaled grid
    xnCVD = seti.T(xnRVD); % seti.T, because we use complex values in this function

    u = G(xnCVD); % u: matrix: nInv^dim
    % u = xnCMD

    % h = seti.hInv
    % N = seti.nInv
    % G = seti.G (reshape vector to matrix)

    % test this code:
    % 1. run start (to set some vectors and functions)
    % 2. gradientNeumann(seti.S(seti.qROIexact),seti.hInv,seti.nInv,seti.G,seti)

Compare slow and fast version of computation of gradient u

    if 0
        disp('test gradient: slow and fast version')
        for dim = 2:3

            fprintf('dimension %g\n',dim)
            N = 100;
            if dim == 2
                u = rand(N,N);
            elseif dim == 3
                u = rand(N,N,N);
            end
            h = 0.1;

            tic
            disp('  gradient u calculation (slow version)')
            gu = guCalcSlow(u,h,N);
            toc

            tic
            disp('  gradient u calculation (fast version)')
            gu2 = guCalc(u,h,N);
            toc

            fprintf('  max difference: %g \n',max(unique(abs(gu-gu2))))
        end
        error('stop: test in gradient neumann end')
    end

Computation of the gradient

    gu = guCalc(u,h,N); % gu is a COMPLEX array of size
                        % 2 x N x N in case of 2D
                        % and size 3 x N x N x N in case of 3D.

Save as real array

    res = seti.S(gu);   % identify complex C with real R x R, then res is real
                        % of size 2*2 x N x N in case of 2D
                        % and size 2*3 x N x N x N in case of 3D.

    end

Code: subfunction: guCalc (fast computation)

    function res = guCalc(u,h,N)
    dim = length(size(u));

    if dim == 2
        gu = zeros(dim,N,N);

        % gu(1,:,:) = grad(u)_{i,j}^1 = 1/h*(u(i+1,j)-u(i,j)) if i < N (0 else)
        uRowShift = circshift(u,[-1 0]);
        uRowShift(end,:) = uRowShift(end-1,:); % to get 0 entry if i = N
        gu(1,:,:) = 1/h*(uRowShift-u);

        % gu(2,:,:) = grad(u)_{i,j}^2 = 1/h*(u(i,j+1)-u(i,j)) if j < N (0 else)
        uColShift = circshift(u,[0 -1]);
        uColShift(:,end) = uColShift(:,end-1); % to get 0 entry if j = N
        gu(2,:,:) = 1/h*(uColShift-u);

        res = gu;
    elseif dim == 3
        gu = zeros(dim,N,N,N);

        % gu(1,:,:,:) = grad(u)_{i,j,k}^1 = 1/h*(u(i+1,j,k)-u(i,j,k)) if i < N (0 else)
        u1Shift = circshift(u,[-1 0 0]);
        u1Shift(end,:,:) = u1Shift(end-1,:,:); % to get 0 entry if i = N
        gu(1,:,:,:) = 1/h*(u1Shift-u);

        % gu(2,:,:) = grad(u)_{i,j,k}^2 = 1/h*(u(i,j+1,k)-u(i,j,k)) if j < N (0 else)
        u2Shift = circshift(u,[0 -1 0]);
        u2Shift(:,end,:) = u2Shift(:,end-1,:); % to get 0 entry if j = N
        gu(2,:,:,:) = 1/h*(u2Shift-u);

        % gu(3,:,:) = grad(u)_{i,j,k}^3 = 1/h*(u(i,j,k+1)-u(i,j,k)) if k < N (0 else)
        u3Shift = circshift(u,[0 0 -1]);
        u3Shift(:,:,end) = u3Shift(:,:,end-1); % to get 0 entry if k = N
        gu(3,:,:,:) = 1/h*(u3Shift-u);

        res = gu;
    else
        error('gradientNeumann.m: dimension must be 2D or 3D.')
    end
    end

Code: subfunction: guCalcSlow (slow computation)

Is not used any more.

    function res = guCalcSlow(u,h,N)

    dim = length(size(u));
    if dim == 2
        gu = zeros(dim,N,N); % grad(u) stored in components gu(1,:,:), gu(2,:,:), (3D: gu(3,:,:))
        % 2D
        for l = 1:dim
            gu(l,:,:) = gradulCalc2D(l,u,h,N); % (grad u)_{i,j}^l
        end
        res = gu;
    elseif dim == 3
        gu = zeros(dim,N,N,N);
        % 3D
        for l = 1:dim
            gu(l,:,:,:) = gradulCalc3D(l,u,h,N); % (grad u)_{i,j,k}^l
        end
        res = gu;
    else
        error('gradientNeumann.m: dimension of u must be 2D or 3D.')
    end
    end

Code: subfunction: gradulCalc2D

gradulCalc2D is required in guCalcSlow.

    function res = gradulCalc2D(l,u,h,N)
    % compute (grad u)_{i,j}^l
    % l = 1 and 2 in 2D
    % l = 1, 2, 3 in 3D
    gradul = zeros(N,N);
    for i = 1:N
        for j = 1:N
            if (l == 1 && i < N) || (l == 2 && j < N)
                switch l
                    case 1
                        gradul(i,j) = (u(i+1,j)-u(i,j))/h; % first component
                    case 2
                        gradul(i,j) = (u(i,j+1)-u(i,j))/h; % second component
                end
            else %case i = N
                gradul(i,j) = 0;
            end
        end
    end
    res = gradul;
    end

Code: subfunction: gradulCalc3D

gradulCalc3D is required in guCalcSlow

    function res = gradulCalc3D(l,u,h,N)
    % compute (grad u)_{i,j}^l
    % l = 1 and 2 in 2D
    % l = 1, 2, 3 in 3D
    gradul = zeros(N,N,N);
    for i = 1:N
        for j = 1:N
            for k = 1:N
                if (l == 1 && i < N) || (l == 2 && j < N) || (l == 3 && k < N)
                    switch l
                        case 1
                            gradul(i,j,k) = (u(i+1,j,k)-u(i,j,k))/h; % first component
                        case 2
                            gradul(i,j,k) = (u(i,j+1,k)-u(i,j,k))/h; % second component
                        case 3
                            gradul(i,j,k) = (u(i,j,k+1)-u(i,j,k))/h; % second component
                    end
                else %case i = N
                    gradul(i,j,k) = 0;
                end
            end
        end
    end
    res = gradul;
    end