Login [Register]
Don't have an account? Register now to chat, post, use our tools, and much more.
On a whim, I wrote a program in matlab to simulate the modes of a 2D circular membrane vibrating in 3D space, also known as a drum. The program takes the two variables for the mode of the membrane (n,m), where n is the number of radial nodes, and m is the number of concentric nodes (1=none only center and edge, 2=center, edge, and 1 more, etc), then generates a plot of the membrane, using Laplacian smoothing to create the surface. However, it still has a significant number of flaws, so if anyone would like to mess with my code and see if they can improve it a bit, I'd appreciate the help. Check out some sample plots and the code itself below.

drumsolve(8,8);


drumsolve(2,6);


drumsolve(6,2);


drumsolve(2,1);



Code:
% Christopher Mitchell, www.cemetech.net
% Circular Membrane (Drum) Mode Solver/Simulator

function [] = drumsolve(n,m)

gran_rad=0.02;
gran_theta=(2*pi)/360;
changetol = 0.0001;

radmax = 5;
zmax = 1;

rad=[0:gran_rad:radmax];
theta=[0:gran_theta:2*pi];

%generate the polar mesh
for a=1:length(rad)
    x(:,a)=rad(a)*cos(theta);
    y(:,a)=rad(a)*sin(theta);
    movable(1:length(theta),a)=1;
end

%set up immovable radii
if (n>0)
    for a=1:n
        myangle=(a-1)*(pi./n);
        movable(1+round(myangle./gran_theta),:)=0;
        z(1+round(myangle./gran_theta),:)=0;
        movable(1+round((pi+myangle)./gran_theta),:)=0;
        z(1+round((pi+myangle)./gran_theta),:)=0;
    end
end

% set up immovable concentric circles
if (m>0)
    for a=1:m+1
        myrad=(a-1)*(radmax/m);
        movable(:,1+round(myrad./gran_rad))=0;
        z(:,1+round(myrad./gran_rad))=0;
    end
end

%set up initial z values to smooth
rowstart = 1;

for b=2:length(rad)-1
    if (m==0 || movable(2,b)==1)
        blockstart = rowstart;
        for a=1:length(theta)
            if movable(a,b)==1
                z(a,b)=blockstart;
            else
                blockstart=blockstart*-1;
            end
        end
    else
        rowstart = rowstart*-1;
    end
end

%smooth it out using Laplace's method
done=0;
while done==0
    done=1;
    for b=1:length(rad)
        for a=1:length(theta)
            if (movable(a,b) == 1)
                div = 3;
                if (a==1)
                    neighbortotal = z(length(theta),b);
                else
                    neighbortotal = z(a-1,b);
                end
                if (a==length(theta))
                    neighbortotal = neighbortotal+z(1,b);
                else
                    neighbortotal = neighbortotal+z(a+1,b);
                end
                if (b>1)
                    neighbortotal = neighbortotal+z(a,b-1);
                else
                    neighbortotal = neighbortotal+z(mod((a+round(pi./gran_theta)-1),length(theta))+1,2);
                end
                if (b<length(rad))
                    neighbortotal = neighbortotal+z(a,b+1);
                    div = div+1;
                end
                new = neighbortotal./div;
                if abs(new-z(a,b))>changetol
                    done = 0;
                end
                z(a,b)=new;
            end
        end
    end
end
z=z.*zmax./max(max(abs(z)));
figure(1);
for i=1:5;
    for rho=0:pi./12:2*pi;
        coeff = sin(rho);
        clf;
        mesh(x,y,coeff.*z);
        axis([-radmax radmax -radmax radmax -2*zmax 2*zmax -1*zmax zmax]);
        M(1)=getframe();
    end
end
oh snap, I've always wanted one of these Wink

seriously though, I'd be happy to test it if
a) I knew enough about circular membranes with concentric and radial nodes
b) I had a matlab license.
No worries. It currently fails at two things:
1) there's a discontinuity at the boundaries between 0 and 360 degrees (it's rectangular mapped to polar)
2) It doubles the mode values for m because I misread my notes from Musical Instrument Design class. No biggie.
Which version of Matlab are you running? I had to comment out the line with "clf;" in order for the movie to work. I'm using v7.6.0 (R2008a).

Also, I believe the method of the "smooth function" is very similar to a numerical laplacian, but skewed and streched to polar coordinates. Check out:
http://en.wikipedia.org/wiki/Hearing_the_shape_of_a_drum
zeromp wrote:
Which version of Matlab are you running? I had to comment out the line with "clf;" in order for the movie to work. I'm using v7.6.0 (R2008a).

Also, I believe the method of the "smooth function" is very similar to a numerical laplacian, but skewed and streched to polar coordinates. Check out:
http://en.wikipedia.org/wiki/Hearing_the_shape_of_a_drum
Welcome to Cemetech. Yes, I did indeed perform a Numerical Laplacian in the discretized domain, but I figured no one here would know what I was talking about if I said that. I'm running 2007a at the moment; I'm surprised the movie didn't work.
It's been a while since I've done my PDEs, but I think that may be useful here...although you did an awesome job already.
zeromp wrote:
It's been a while since I've done my PDEs, but I think that may be useful here...although you did an awesome job already.
Oh thanks, I plan to fix up the bugs and repost soon, but I'd definitely appreciate the help if you have a chance to look into it. Smile May I ask what brought you to this?
I was on ticalc.org and noticed the TI-83 programming contest so I entered that and ever since, I've been intrigued with the projects that you work on. I'm guessing you're an EE, but I'm not 100% sure.

Anyway, once I saw that you like to mess around with Matlab, I had to check it out. I do quite a bit of stuff on Matlab myself and that's what brought me to this.
zeromp wrote:
I was on ticalc.org and noticed the TI-83 programming contest so I entered that and ever since, I've been intrigued with the projects that you work on. I'm guessing you're an EE, but I'm not 100% sure.

Anyway, once I saw that you like to mess around with Matlab, I had to check it out. I do quite a bit of stuff on Matlab myself and that's what brought me to this.
That's correct, I'm a senior undergraduate EE. I'm glad to hear you are intrigued by my projects, and it would be neat if you'd post some of your own Matlab projects
*bump* Added some fixes and scaling factors; I notice it causes some very weird behavior, not sure if that's in error or actually an accurate simulation.


Code:
% Christopher Mitchell EE 2009
% Circular Membrane Mode Solver/Simulator v1.1

function [] = drumsolve(n,m)

figure(1);

gran_rad=0.1;
gran_theta=(2*pi)/360;
changetol = 0.001;

radmax = 5;
zmax = 1;

rad=[0:gran_rad:radmax];
theta=[0:gran_theta:2*pi];

%generate the polar mesh
for a=1:length(rad)
    x(:,a)=rad(a)*cos(theta);
    y(:,a)=rad(a)*sin(theta);
    movable(1:length(theta),a)=1;
end

%set up immovable radii
if (n>0)
    for a=1:n
        myangle=(a-1)*(pi./n);
        movable(1+round(myangle./gran_theta),:)=0;
        z(1+round(myangle./gran_theta),:)=0;
        movable(1+round((pi+myangle)./gran_theta),:)=0;
        z(1+round((pi+myangle)./gran_theta),:)=0;
    end
end

% set up immovable concentric circles
if (m>0)
    for a=1:m-1
        myrad=radmax-(a-1)*(radmax/(m-1));
        movable(:,1+round(myrad./gran_rad))=0;
        z(:,1+round(myrad./gran_rad))=0;
    end
end

%set up initial z values to smooth
rowstart = 1;

for b=2:length(rad)-1
    if (m==0 || movable(2,b)==1)
        blockstart = rowstart;
        for a=1:length(theta)
            if movable(a,b)==1
                z(a,b)=blockstart;
            else
                blockstart=blockstart*-1;
            end
        end
    else
        rowstart = rowstart*-1;
    end
end

%smooth it out using the numerical Laplacian method, modified for polar
%coordinates, or if you prefer, the polar-modified numerical Laplacian
done=0;
iter=0;
while done==0

%    disp(sprintf('min max %f %f',min(min(z)),max(max(z))));

    done=1;
    for b=1:length(rad)
        for a=1:length(theta)
            if (movable(a,b) == 1)
                div = 1;
                neighbortotal = 0;
                if (b ~= 1)
                    rscale = sin((pi-gran_theta)/2)./((b-1).*sin(gran_theta));
                    div = div + 2.*rscale;
                    if (a==1)
                        neighbortotal = rscale.*z(length(theta),b);
                    else
                        neighbortotal = rscale.*z(a-1,b);
                    end
                    if (a==length(theta))
                        neighbortotal = neighbortotal + rscale.*z(1,b);
                    else
                        neighbortotal = neighbortotal + rscale.*z(a+1,b);
                    end
                end
                if (b > 1)
                    neighbortotal = neighbortotal + z(a,b-1);
                else
                    neighbortotal = neighbortotal + z(mod((a+round(pi./gran_theta)-1),length(theta))+1,2);
                end
                if (b<length(rad))
                    neighbortotal = neighbortotal + z(a,b+1);
                    div = div + 1;
                end
                new = neighbortotal./div;
                if abs(new-z(a,b))>changetol
                    done = 0;
                end
                z(a,b)=new;
            end
        end
    end
iter = iter+1;

end

disp(sprintf('Converged in %.0f iterations.',iter));
z=z.*zmax./max(max(abs(z)));
disp(sprintf('min %f max %f',min(min(z)),max(max(z))));

for i=1:5;
    for rho=0:pi./12:2*pi;
        coeff = sin(rho);
        clf;
        mesh(x,y,coeff.*z);
        axis([-radmax radmax -radmax radmax -2*zmax 2*zmax -1*zmax zmax]);
        M(1)=getframe();
    end
end
Whoa, you saved my life with this coding... Seriously, it helps me save time figuring out the whole thing for my End of the Year project for College.

Never used MATLAB yet and I just started using it, but can I use the coding for reference or something? I'm new to the whole open source, code sharing, thing on the net.

If not, can I base my programming with what you provided? Did you discretize Laplacian or did you use the Bessel function for this?
Haha, glad it could help, and welcome to Cemetech. As long as you credit me (Christopher Mitchell) with the original coding, you're welcome to use and modify my code as you see fit, as long as it's not for commercial profit. Smile I used a discretized numerical Laplacian, aka the Monte Carlo iterative solution scheme, so it's an approximated solution. I'm not sure if it's perfectly accurate, but I hope I got most of it down properly. The polar coordinate scheme is slightly kludgy, but it's as good as anything else I could find to do the same thing. What's your application?
Alright then! Yeah, I know when to give credits where it's due, plus since I'm in college in montreal we don't normally see that in our year. In University yea, but not in our years of college. So, if I presented something like this, I don't think I can get through without ending with a fail due to plagirism.

Anyways, what do you mean by application? Thanks again for the work though =D
FuyuKaze wrote:
Alright then! Yeah, I know when to give credits where it's due, plus since I'm in college in montreal we don't normally see that in our year. In University yea, but not in our years of college. So, if I presented something like this, I don't think I can get through without ending with a fail due to plagirism.

I'm assuming you're using the European naming scheme where college == high school, and university == college, or are you making the distinction that college == undergraduate work and university == graduate work?

Quote:
Anyways, what do you mean by application? Thanks again for the work though =D

He means what specific task would you be using it for?
really useful......thanks
miketyson986 wrote:
really useful......thanks
Cheers, welcome to Cemetech. May I ask what it's useful to you for?
  
Register to Join the Conversation
Have your own thoughts to add to this or any other topic? Want to ask a question, offer a suggestion, share your own programs and projects, upload a file to the file archives, get help with calculator and computer programming, or simply chat with like-minded coders and tech and calculator enthusiasts via the site-wide AJAX SAX widget? Registration for a free Cemetech account only takes a minute.

» Go to Registration page
Page 1 of 1
» All times are GMT - 5 Hours
 
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum