Thread Subject:

How do I create the 3D central difference matrices (using kron)?

Say I have a 3x3x3 matrix (tensor?):
a=[1:27]';
A=reshape(a, [3 3 3]);

How do I create the 3D central difference matrices: Dx3, Dy3, Dz3,
such that: reshape(Dx3*a, [3 3 3]): partial derivative wrt x and so on?

I recently learned how to do this in 2D (3x3 matrix):
E = sparse(2:3, 1:2, 1, 3, 3);
D = (E' - E) / 2;
I = speye(3);
Dx = kron(I,D);
Dy = kron(D,I);

But how do I extend this to 3D?

I would prefer a solution using kron and my matrix D from above.

Thanks in advance for any answers!

On Friday, June 29, 2012 11:34:45 AM UTC+2, Andy wrote:
> Say I have a 3x3x3 matrix (tensor?):
> a=[1:27]';
> A=reshape(a, [3 3 3]);
>
> How do I create the 3D central difference matrices: Dx3, Dy3, Dz3,
> such that: reshape(Dx3*a, [3 3 3]): partial derivative wrt x and so on?
>
> I recently learned how to do this in 2D (3x3 matrix):
> E = sparse(2:3, 1:2, 1, 3, 3);
> D = (E' - E) / 2;
> I = speye(3);
> Dx = kron(I,D);
> Dy = kron(D,I);
>
> But how do I extend this to 3D?
>
> I would prefer a solution using kron and my matrix D from above.
>
> Thanks in advance for any answers!

Sorry for the duplicate. Waited 45 min and my first post did not show up.
Andreas

Andy <andreas.paulsen@geocap.no> wrote in message <2ec64a27-81e3-4e6d-aaf3-c2fee75f331c@googlegroups.com>...
> Say I have a 3x3x3 matrix (tensor?):
> a=[1:27]';
> A=reshape(a, [3 3 3]);
>
> How do I create the 3D central difference matrices: Dx3, Dy3, Dz3,
> such that: reshape(Dx3*a, [3 3 3]): partial derivative wrt x and so on?
>
> I recently learned how to do this in 2D (3x3 matrix):
> E = sparse(2:3, 1:2, 1, 3, 3);
> D = (E' - E) / 2;
> I = speye(3);
> Dx = kron(I,D);
> Dy = kron(D,I);
>
> But how do I extend this to 3D?
===============

As below,

E = sparse(2:3, 1:2, 1, 3, 3);
D = (E' - E) / 2;
I = speye(3);
Dx = kron(I,kron(I,D));
Dy = kron(I,kron(D,I));
Dz = kron(D,kron(I,I));

However, if you're eventually going to be working with data sizes much larger than 3x3x3, kron can overflow memory limits, depending on you platform, even though Dx,Dy,Dz are in sparse form. A better option is then to use this class:

http://www.mathworks.com/matlabcentral/fileexchange/25969-efficient-object-oriented-kronecker-product-manipulation

For example, on my 32-bit laptop, the kron approach will not work on the data below and only the KronProd class is viable for me.

N=300;
A=rand(N,N,N);
I=speye(N);

D=interpMatrix([-1 0 1]/2, 2, size(A,1), 1).';

Dx=KronProd({D,1},[1,2,2],size(A));


tic;
 Dx*A;
toc;%Elapsed time is 0.723152 seconds.

 
Note also that I used interpMatrix to compute the 1D differencing matrix D, which is available here

http://www.mathworks.com/matlabcentral/fileexchange/26292-regular-control-point-interpolation-matrix-with-boundary-conditions

kl. 11:34:45 UTC+2 fredag 29. juni 2012 skrev Andy flgende:
> Say I have a 3x3x3 matrix (tensor?):
> a=[1:27]';
> A=reshape(a, [3 3 3]);
>
> How do I create the 3D central difference matrices: Dx3, Dy3, Dz3,
> such that: reshape(Dx3*a, [3 3 3]): partial derivative wrt x and so on?
>
> I recently learned how to do this in 2D (3x3 matrix):
> E = sparse(2:3, 1:2, 1, 3, 3);
> D = (E' - E) / 2;
> I = speye(3);
> Dx = kron(I,D);
> Dy = kron(D,I);
>
> But how do I extend this to 3D?
>
> I would prefer a solution using kron and my matrix D from above.
>
> Thanks in advance for any answers!

Thanks Matt! That worked great. Especially your KronProd class.
Andreas