3D matrix manipulation problems

I struggle with this generally and have come up against a barrier again.

I'm trying to create vibrational mode functions for a rectangular membrane. The equation is simple enough: -

 W_n = sin(x*pi*i/L_x)*sin(y*pi*i/L_y)

Here's where I'm up to: -

L_x             = 27.4e-3;  %membrane width (m)
L_y             = 27.4e-3;  %membrane height (m)
N_x             = 10;       %no. of eigenfreqs considered in x dimension
N_y             = 10;       %no. of eigenfreqs considered in y dimension
N               = N_x*N_y;  %total no. of eigenfreqs considered

x               = (linspace(0.5*L_x/N_x,L_x - 0.5*L_x/N_x,N_x))'; 
                            %membrane x-dimension mapping points
                            %each at the centre of an element
y               = (linspace(0.5*L_y/N_y,L_y - 0.5*L_y/N_y,N_y));  
                            %membrane y-dimension mapping points
                            %each at the centre of an element
x_ref           = 1:length(x);
y_ref           = 1:length(y);
W_n             = zeros(N_x,N_y,N);
W_n(:,:,01)     = (sin(x(x_ref,:,:)*1*pi/L_x)...
                   *sin(y(:,y_ref,:)*1*pi/L_y));
W_n(:,:,02)     = (sin(x(x_ref,:,:)*1*pi/L_x)...
                   *sin(y(:,y_ref,:)*2*pi/L_y));
W_n(:,:,18)     = (sin(x(x_ref,:,:)*2*pi/L_x)...
                   *sin(y(:,y_ref,:)*8*pi/L_y));
W_n(:,:,46)     = (sin(x(x_ref,:,:)*5*pi/L_x)...
                   *sin(y(:,y_ref,:)*6*pi/L_y));            
W_n(:,:,89)     = (sin(x(x_ref,:,:)*9*pi/L_x)...
                   *sin(y(:,y_ref,:)*9*pi/L_y));            
W_n(:,:,100)    = (sin(x(x_ref,:,:)*10*pi/L_x)...
                   *sin(y(:,y_ref,:)*10*pi/L_y));

The numbers on the left hand side - 01, 02, 18, 46, 89, 100 - actually need to be 1,2,3,...,100, i.e.

n               = 1:100;    %mode number counting from 1 to 100

These refer to the mode numbers - although this is a bit odd - but still not too complicated. Basically if n = 01, then the mode is (1,1). If n = 59 the mode is (6,9), i.e. i = 6 & j = 9. I have solved this using the following: -

i               = floor((n/10 - n/1000)) + 1;       
                            %mode number from 1 to 10 in x dimension
j               = n - 10*(floor((n/10 - n/1000)));  
                            %mode number from 1 to 10 in y dimension

Now I just need to put n, i and j into my W_n equation, but I can't figure out how!

Any help would be greatly appreciated :)

Answer 1

[Simon]

Hi!

I don't understand why the x and y spacing depends on the number of modes/frequencies.

Take a look at the following code:

L_x             = 27.4e-3;  %membrane width (m)
L_y             = 27.4e-3;  %membrane height (m)
N_x             = 5;       %no. of eigenfreqs considered in x dimension
N_y             = 10;       %no. of eigenfreqs considered in y dimension
N               = N_x*N_y;  %total no. of eigenfreqs considered

NumX = 50; % number of mapping points in x-direction
NumY = 50; % number of mapping points in y-direction

% create X and Y array in 2d
[X,Y] = meshgrid(...
    (linspace(0.5*L_y/NumY, L_y - 0.5*L_y/NumY, NumY)), ...
    (linspace(0.5*L_x/NumX, L_x - 0.5*L_x/NumX, NumX)));

% modify X and Y array for 3d
XFull = repmat(X, [1 1 N]);
YFull = repmat(Y, [1 1 N]);

% create mode array for X
Mx = ones(NumX, NumY, N_y);
MxFull = [];
for n = 1:N_x
    MxFull = cat(3, MxFull, n*Mx);
end

% create mode array for Y
My = ones(NumX, NumY);
MyFull = [];
for n = 1:N_y
    MyFull = cat(3, MyFull, n*My);
end
MyFull = repmat(MyFull, [1 1 N_x]);

% modes for each direction
Wx_n = arrayfun(@(x,mx) sin(x .* (mx *pi/L_x)), XFull, MxFull);
Wy_n = arrayfun(@(y,my) sin(y .* (my *pi/L_x)), YFull, MyFull);

% superposition of modes
W_n = Wx_n .* Wy_n;

% plot
figure(1); cla
surface(W_n(:,:,23));

[Tom]

Sorry - yeah I didn't explain that. I don't know for sure (I'm replicating this theory from a few papers on the subject), but I'm pretty sure it is basically that if the maximum no. of modes is (10,10), then the least number of elements needed to give a basic representation of this with any degree of accuracy would be an mesh made up of 10x10 elements. The x-y points are points at the centre of each element.

I'll check out the code now...

[Simon]

I agree with the number of elements vs. number of modes. It's kind of sampling theorem of sin function. But keep in mind that it is the least number. Compare the code above for

N_x = 10;
N_y = 10;
NumX = 10;
NumY = 10;

and

N_x = 10;
N_y = 10;
NumX = 100;
NumY = 100;

Look at

surface(W_n(:,:,100));

[Tom]

That is amazing thank you so much. One thing - when the plot seems to start from 1 in both x and y dimensions rather than 0. Is there a way of getting it to from 0?

[Simon]

Use the real x and y coordinates of the mapping points

figure(1); cla
surface((linspace(0.5*L_y/NumY,L_y - 0.5*L_y/NumY,NumY)), ...
    (linspace(0.5*L_x/NumX,L_x - 0.5*L_x/NumX,NumX)), ...
    W_n(:,:,100));
axis tight

[Tom]

That's perfect - thank you!!

[Tom]

Hi again Simon, I'm just getting into this again and I've noticed that when I examine the 1st mode, it appears off-centre - whereas it should be completely central on the membrane.

This can be seen when the end of the surface command is

 W_n(:,:,1);

Do you know why this is?

Many thanks,

Tom

[Tom]

I think I may have solved this by using

 (linspace(L_y/NumY, L_y, NumY))

instead of

 (linspace(0.5*L_y/NumY,L_y - 0.5*L_y/NumY,NumY))

Do you think this solution is correct?

Thanks