I am getting Error: Error using * Arguments must be 2-D, or at least one argument must be scalar.

I have a 3D function where I am testing taking derivative along x,y, and z direction. My issue is that taking derivative wrt z is giving an error

clearvars; clc; close all;
Nx = 4;
Ny = 4;
Nz = 4;
%-----
Lx = 2*pi; %8; %128;
Ly = 2*pi; 
% Set the number of grid points 

%Set-up grids:
x   = (0:Nx-1)/Nx*2*pi; 
y   = (0:Ny-1)/Ny*2*pi; 

kx = (2*pi/Lx)*([0:Nx/2-1 , -Nx/2:-1]);
ky = (2*pi/Ly)*([0:Ny/2-1 , -Ny/2:-1]);
zgl = -cos(pi*(0:Ny)/Ny)'; %Gauss-Lobatto chebyshev points 
[X,Y,Z] = meshgrid(x,y,zgl);

%Chebyshev matrix for Gauss-Lobatto
VGL = cos(acos(zgl(:))*(0:Nz));
dVGL = diag(1./sqrt(1-zgl.^2))*sin(acos(zgl)*(0:Nz))*diag(0:Nz);
dVGL(1,:)    = (-1).^(1:Nz+1).*(0:Nz).^2;
dVGL(Nz+1,:) = (0:Nz).^2;
%Diferentiation matrix for Gauss-Lobatto points
Dgl = dVGL/VGL;
D   = Dgl; %first-order derivative matrix
%------------------
ubar = Z.^2 .* sin( (2*pi / Lx) * X) .* sin( (2*pi / Ly) * Y) ; %Y^2 * sin( (2*pi / Lx) * x); 
uh = fft(fft(ubar,[],2),[],1); %fft along x and y only

duhdxk = derivk(uh, kx); %works
duhdx = real(ifft(ifft(duhdxk,[],2),[],1));

exactDx = (2*pi)/(Lx)* Z.^2 .* cos( (2*pi / Lx) * X) .* sin( (2*pi / Ly) * Y);
DEx = exactDx -duhdx; 
% y derivative 
duhdyk = derivk(uh, ky'); %works
duhdy = real(ifft(ifft(duhdyk,[],2),[],1));

exactDy = (2*pi)/(Ly)* Z.^2 .* sin( (2*pi / Lx) * X) .* cos( (2*pi / Ly) * Y);
DEy = exactDy - duhdy;
%%
%z derivative 
for e = 1:Nz
    duhdzk(:,:,e) = D * uh(:,:,e); % ERROR: D is 5x5 size and uh is 4x4x5 size
end
exactDz = 2 .* Z .* sin( (2*pi / Lx) * X) .* sin( (2*pi / Ly) * Y);

I am aware the error comes from multiplying D an uh with the incorrect sizes so I tried creating a for loop to change the uh size to 4x5 but it's not working.

Full error:

Error using  * 
Arguments must be 2-D, or at least one argument must be scalar. Use TIMES (.*) for elementwise
multiplication, or use PAGEMTIMES to apply matrix multiplication to the pages of N-D arrays.

Error in FourierCheby3D (line 77)
    duhdzk(:,:,e) = D * uh(:,:,e);
 

[KSSV]

D * uh(:,:,e); % ERROR: D is 5x5 size and uh is 4x4x5 size

How you can multiply 5x5 matrix and 4x4 matrix? Go through the basics of matrix multiplication.

[Walter Roberson]

We do not have derivk() to test the code with.

[Jamie Al]

Yes, sorry this is the function:

function data_deriv = derivk(fk,k)

% Takes derivative of a 2D matrix using Fourier transforms

data_deriv  = 1i * k .* fk;

[Jamie Al]

I know I can't multiply 5x5 matrix by 4x4, that's why I was trying somehow to change the indexing to 4x5 maybe. I am not sure how to take the derivative wrt z otherwise.

[Walter Roberson]

I do not get the reported message. I get

Error using *

Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform elementwise multiplication,

use '.*'.

Error in q1690600 (line 45)

duhdzk(:,:,e) = D * uh(:,:,e); % ERROR: D is 5x5 size and uh is 4x4x5 size

You could index to 4 x 5 by using

 duhdzk(:,:,e) = D(:,1:4) * uh(:,:,e);

or

  duhdzk(:,:,e) = D(:,2:5) * uh(:,:,e);

or

   duhdzk(:,:,e) = D(:,setdiff(1:5,randi(5))) * uh(:,:,e);

As outside observers, we have no information as to what what make the most sense in the situation.

VGL = cos(acos(zgl(:))*(0:Nz));
dVGL = diag(1./sqrt(1-zgl.^2))*sin(acos(zgl)*(0:Nz))*diag(0:Nz);
dVGL(1,:)    = (-1).^(1:Nz+1).*(0:Nz).^2;
dVGL(Nz+1,:) = (0:Nz).^2;

Maybe those should be 0:Nz-1 and 1:Nz ?

[Jamie Al]

First of all thanks for the help.

I don't know if you're able to run my code, but you can see I am trying to compare numerical derivatives to exact ones. So, my x and y derivatives are matching (taking derivatives in Fourier space). But, the third derivative along z is creating an issue. I am taking the derivative along z using chebyshev derivative matrix D which usually has a size of Nz+1 x Nz+1. While, your suggestions work, now I can't compare between my exact derivative and the numerical one. So, I get the error:

Array dimensions must match for binary array op.

Error in FourierCheby3D (line 86)
DEz = exactDz - duhdz;

I guess I am not appraoching this correctly.

This comes from the definition of the chebyshev matrix and I am following with a textbook, so I am not sure if I should change things here:

VGL = cos(acos(zgl(:))*(0:Nz));
dVGL = diag(1./sqrt(1-zgl.^2))*sin(acos(zgl)*(0:Nz))*diag(0:Nz);
dVGL(1,:)    = (-1).^(1:Nz+1).*(0:Nz).^2;
dVGL(Nz+1,:) = (0:Nz).^2;

[Walter Roberson]

https://scicomp.stackexchange.com/questions/24470/how-do-i-form-the-chebyshev-differentiation-matrix-in-matlab

Answer 1

[Jamie Al]

So I was able to resolve this issue, by indexing in a very specific way:

[X,Z,Y] = meshgrid(x,zgl,y); %now this is a Z-by-X-by-Y 3D grid
%Taking FFT of X and Y is changed to 2nd and 3rd dimension instead
uh = fft(fft(ubar,[],2),[],3);
%Thus, x derivative is
duhdxk = derivk(uh, kx);
%y derivative
duhdyk = derivk(uh, reshape(ky, [1,1,Ny])); 
%z derivative 

for e = 1:Nz 
    duhdzk(:,:,e) = D * uh(:,:,e);
end
duhdz = real(ifft(ifft(duhdzk,[],2),[],3));

This works, however, I am not sure if I can index things correctly with a

[X,Y,Z] = meshgrid(x,y,zgl);