4 DOF Model

[Steve]

Jim, I think both John and Dom have been through the derivations (they will correct me if I'm wrong!) and I think they got what's in the book. John has also managed to get pretty similar charts to mine, here: viewtopic.php?f=33&t=4636&start=50#p55463, so at this point I'd say have another look and see if there's an error in your workings. A Matlab solution will certainly be good to have, as Excel is pretty flaky on this stuff, so I'm looking forward to seeing what you have, when you're ready.

I was working through this stuff yesterday too (up to the point of the y's - haven't moved to the next phase yet). I used Excel (both 2003 & 2007 - both got the same answers), and came up with exactly the same solutions (including Dbar) that John reported in the screenshot at the bottom of his post here: viewtopic.php?f=33&t=4636&start=25#p54820. The only big difference that I've spotted between his excel & VB solution was the ya value. I've re-worked that calc about a dozen times & keep arriving at the same answer (excel version). For what it's worth, all of my calculations have been single, nested functions - I was getting the same answers when they were split up.
Somehow we're all misinterpreting your equations (or making the same mistakes) to come up with the same, incorrect answers. The bit that's really doing my head in is that John's final VB version came up with radically different values to your 100Hz intermediates (viewtopic.php?f=33&t=4636&start=25#p54844), but then arrive at the same/similar graph! Do you have any thoughts that you could offer on this? I haven't used VB in a few years, but perhaps it's the only way forward?
Thanks,
Steve

[Trevor Gore]

Somehow we're all misinterpreting your equations (or making the same mistakes) to come up with the same, incorrect answers. The bit that's really doing my head in is that John's final VB version came up with radically different values to your 100Hz intermediates (viewtopic.php?f=33&t=4636&start=25#p54844), but then arrive at the same/similar graph! Do you have any thoughts that you could offer on this?

I have nothing more to offer at the moment, sorry. All I can tell you is that the Excel solutions are pretty flaky and that I'm hanging out for Jim's Matlab solution, as that is likely to be the most robust to date and should cast a lot of light on where/how our solutions are diverging.

[Steve]

I have nothing more to offer at the moment, sorry. All I can tell you is that the Excel solutions are pretty flaky and that I'm hanging out for Jim's Matlab solution, as that is likely to be the most robust to date and should cast a lot of light on where/how our solutions are diverging.

Bummer. I'll keep plugging away then.
Thanks,
Steve

[Jim Kirby]

OK, here is what I get. First, two pictures. The first is the book example with the sign change I claimed above needs to be made. The second is what I get without the sign change (i.e., using the formula as given in the book.)

"correct?" version

fig1.jpg (30.64 KiB) Viewed 53612 times

original version

fig2.jpg (30.17 KiB) Viewed 53612 times

[Jim Kirby]

Here is the matlab .m file. Please let me know if you find any errors.

Jim

[Jim Kirby]

I couldn't find a file extension that the board would let me upload (what do you use for ASCII text files?
Here is the code inline.

Jim


%--------------------------------------------------------------------------
%
% gore_gillet_4dof.m
%
% Implements the 4 DOF model described in section 2.4.1 in Gore and Gillet,
% etc.
%
%
% copyright Jim Kirby, 11/27/12 - 12/28/12, free for use by anyone.
%
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%
% Physical constants(Table 2.4-1, pg 2-37.
%
%--------------------------------------------------------------------------

mt=0.07; % mass of top
mb=0.094; % mass of back
ms=0.8; % mass of sides
ma=0.000591; % sound hole air mass

At=0.0430; % area of top
Ab=0.04; % area of back
Aa=0.0063; % soundhole area
As=0.025; % area of sides

Rt=3.0; % damping coefficient, top
Rb=7.0; % damping coefficient, back
Rs=10.0; % damping coefficient, sides
Ra=0.014; % damping coefficient, air

Kt=68000; % equivalent top plate stiffness
Kb=130000; % equivalent back plate stiffness

c2=340^2; % speed of sound squared
rho=1.205; % density of air
V=0.015; % volume of box

pref=50;

kappa=c2*rho/V;

%--------------------------------------------------------------------------
%
% Uncoupled natural frequencies of the four modes
%
%--------------------------------------------------------------------------

omt=sqrt( (Kt + kappa*At^2)/mt);

oms=sqrt( (Kb + Kt)/ms);

omb=sqrt( (Kb + kappa*Ab^2)/mb);

oma=sqrt( kappa*Aa^2/ma);

%--------------------------------------------------------------------------
%
% Coupling constants
%
%--------------------------------------------------------------------------

alphat=kappa*Aa*At;

alphab=kappa*Aa*Ab;

alphbt=kappa*Ab*At;

%--------------------------------------------------------------------------
%
% frequency, angular frequency
%
%--------------------------------------------------------------------------

f=(0:0.1:500);

om=2*pi*f;

%--------------------------------------------------------------------------
%
% Coefficients in final solution
%
%--------------------------------------------------------------------------

Dt=mt*(omt^2 - om.^2) + i*Rt*om;

Ds=ms*(oms^2 - om.^2) + i*Rs*om;

Db=mb*(omb^2 - om.^2) + i*Rb*om;

Da=ma*(oma^2 - om.^2) + i*Ra*om;

Dbar= Dt.*Ds.*Db.*Da - 2*alphbt*Kt*Kb*Da + 2*alphbt*alphab*alphat*Ds + ...
2*alphat*alphab*Kt*Kb -(alphab^2)*(Dt.*Ds - (Kt^2)) ...
-(alphat^2)*(Ds.*Db - (Kb^2)) - (alphbt^2)*Ds.*Da - (Kb^2)*Dt.*Da ...
-(Kt^2)*Db.*Da;

%--------------------------------------------------------------------------
%
% final displacements
%
%--------------------------------------------------------------------------

F=0.4; % external applied force amplitude

yt=(F./Dbar).*(Ds.*Db.*Da - (alphab^2)*Ds - (Kb^2)*Da);

ys=(F./Dbar).*(Kt*Db.*Da - Kt*(alphab^2) + alphbt*Kb*Da - alphat*alphab*Kb);

yb=(F./Dbar).*(-Kt*Kb*Da - alphbt*Da.*Ds + alphab*alphat*Ds);

ya=(F./Dbar).*(Kt*Kb*alphab + alphbt*alphab*Ds - alphat*Ds.*Db + alphat*(Kb^2) );
% The book has a minus sign in front of the first term

%--------------------------------------------------------------------------
%
% Resulting sound pressure (Equ. 2.4-11,12)
%
%--------------------------------------------------------------------------

r=1;

p=(om.*om*rho/(4*pi*r)).*(At*yt + As*ys + Ab*yb + Aa*ya);

db=20*log10(abs(p)/pref);

%--------------------------------------------------------------------------
%
% Plot results
%
%--------------------------------------------------------------------------

plot(f,db), grid on, xlabel('f (hz)'), ylabel('sound pressure level (dB)'),axis([0 500 -100 -20])

[GregL]

Hi Jim,

No answers, I'm afraid. I haven't used matlab before, nor have I tried to solve any of the "book" equations. However, if you could clear up a couple of points, it would improve my understanding and maybe help to find a solution.

- what does an ellipsis mean in your code? Is it something you have edited into the text to make it more readable or does matlab use it?

- what is the "dot notation"? eg:

yt=(F./Dbar).*(Ds.*Db.*Da - (alphab^2)*Ds - (Kb^2)*Da);

what does "F." mean?

- suggestion 1) - put parentheses around absolutely everything. I know that these sorts of programs have an official "order of operations", but I have found through bitter experience that the order of operations may not be strictly respected deep down in the bowels of the program - using parentheses forces correct order of operations;

- suggestion 2) - break big formulae down into small chunks. eg instead of:

yt=(F./Dbar).*(Ds.*Db.*Da - (alphab^2)*Ds - (Kb^2)*Da);

maybe you could write:

yta = (F./Dbar).;
ytb = Ds.*Db.*DDa;
ytc = (alpha^2)*Ds;
ytd = (Kb^2)*Da;
yt = yta * (ytb - ytc - ytd);

Differences may be easier to spot this way.

Thanks,
GregL.

[Jim Kirby]

Greg,

The code got a little mangled by the upload as inline text - the forum software deleted all blanks inside comments or after major punctuations. I'd be glad to send a copy of the original to anyone who sends me an email address. The inline version would run but it's harder to read.

Matlab uses ellipsis to indicate line continuations.

I'm not too worried about order of operation - Matlab is pretty robust that way. I actually tend to overuse parentheses for readability purposes, much beyond what the code actually needs. I'm pretty experienced with it, and am in the process of building an entire toolbox to accompany a text on ocean time series analysis that I'm writing.

Same for breaking down steps. At this point I'm really not too concerned about errors in my own code (but, of course, stuff happens) - I'm more concerned about the sign change w/r G&G's printed equation, which I suspect is a typo - I don't think you can get the figures they show with the sign as it is in their equation, and the step from appendix equation 38 to 43 is very simple and would indicate that 43 has an error. I've been hoping that Trevor would weigh in about this. I re-did the entire derivation myself and get exactly the same results up to the solution forms in terms of determinants, and if 38 is correct, then 43 is wrong.

For those who have an interest in Matlab but without the corresponding required bulge in the wallet (which is steep - boy do these guys have a great marketing strategy - give it to all university people almost for free and then watch them suffer withdrawal when they graduate or leave the hallowed halls), give Octave a try. It's a free Matlab clone that is usually only a version or two behind. Simple codes like mine here would run in Octave with no problem. It only gets hairy when you have complex graphical constructs or special functions from add-on toolboxes.

Jim

[Jim Kirby]

Greg,

I forgot to add an explanation of the .(operator) construct. This indicates an element by element operation on an array where you are retaining the identity of each element. For example, the F./Dbar, where F is a scalar and Dbar is a vector, gives a vector with elements

F/Dbar(1), F/Dbar(2), F/Dbar(3), etc

with the same length as Dbar. It's like looping over the array indices and performing the scalar operation on each element of the array.

Jim

[Jim Kirby]

To summarize, in hopes that someone is checking the math too:

It's not the programming, it's that, if Equ. AII 2-38 is correct (which I think it is), then Equ. AII 2-43 and it's corresponding copy in the main text is wrong. And I suspect the error is purely a typo in the text rather than something that has been incorporated in Trevor's code, which I suspect is correct as it gets the proper form for the main air peak. (Trevor please check your code).

I feel bad causing such a ruckus within my first 10 posts on the forum.

Jim

[Steve]

To summarize, in hopes that someone is checking the math too:

It's not the programming, it's that, if Equ. AII 2-38 is correct (which I think it is), then Equ. AII 2-43 and it's corresponding copy in the main text is wrong. And I suspect the error is purely a typo in the text rather than something that has been incorporated in Trevor's code, which I suspect is correct as it gets the proper form for the main air peak. (Trevor please check your code).

I feel bad causing such a ruckus within my first 10 posts on the forum.

Jim

Jim, I'm a bit rusty on matrices so haven't been able to check the derivations that you have been looking at. What I have done is complete my workings in Excel. My results appear to be the same as your matlab outputs (including your modified AII 2-38)...

Orig.jpg (126.28 KiB) Viewed 53598 times

...& your modified version of AII 2-43

Jim'sVersion.jpg (124.36 KiB) Viewed 53598 times

Like I said, I haven't checked the equation derivations, but your change certainly produces a closer match to Fig 2.4-4, apart from the differences in the 0-100 Hz range.
Hope this helps,
Steve

[Jim Kirby]

Steve - exactly the same as mine.

As for the math, the Wikipedia article on determinants is enough to go by. (I admittedly checked myself).

As for the change in amplitude of the peak, a lot of the peaks near resonance are based on dividing small numbers by small numbers, and can be really sensitive to roundoff errors.

[Steve]

As for the math, the Wikipedia article on determinants is enough to go by. (I admittedly checked myself).

I had a go at this & have to agree with you that I'm getting a positive value for the first term in in Equ AII 2-43 (ie., F(+Kt.Kb.alphaab +...). My solutions for AII 2-40 & 2-41 were the same as documented in the book, but I'm getting a different answer for AII 2-42 - for the final term on the top line, I'm getting Ds.alphaab^2, rather than Ds.alphaab.alphaat. Given my skill level in the area of matrices I'm happy to be corrected. For what it's worth, I changed my calcs for yb to see if it made a difference, but alas the numbers are similar enough that the output graphs looked exactly the same, so right or wrong, it's not going to make a difference.
Steve

[Trevor Gore]

This looks like real progress, so thanks everyone.

I feel bad causing such a ruckus within my first 10 posts on the forum.

No problem at all, Jim. This is what this thread is all about.

I won't have time over the next couple of days to get into this in detail, but please don't stop posting just because I'm not saying anything. If John and Dom have any further insights, please chip in, too. I'll go through my derivations and code again and see what I find and will get back.

Jim, thanks for the tip about Octave.

[Jim Kirby]

As for the math, the Wikipedia article on determinants is enough to go by. (I admittedly checked myself).

I had a go at this & have to agree with you that I'm getting a positive value for the first term in in Equ AII 2-43 (ie., F(+Kt.Kb.alphaab +...). My solutions for AII 2-40 & 2-41 were the same as documented in the book, but I'm getting a different answer for AII 2-42 - for the final term on the top line, I'm getting Ds.alphaab^2, rather than Ds.alphaab.alphaat. Given my skill level in the area of matrices I'm happy to be corrected. For what it's worth, I changed my calcs for yb to see if it made a difference, but alas the numbers are similar enough that the output graphs looked exactly the same, so right or wrong, it's not going to make a difference.
Steve

Steve,

Ah-ha, you found another typo in AII 2-37. That lower left element in the 4 by 4 should always be alpha_at except in y_t, where that is the column that is replaced. With a proper alpha_at in AII 2-37, AII 2-42 is correct.

[Dominic]

Hey guys, glad the batton has been picked up. I have been approaching the problem from the angle of defeating excels stupid treatment of small numbers. I have also discovered that the Polish version of excel is worse by several decimal places.
Anyway, splitting up the equations into smaller parts before resum would theoretically solve the problem. Got it half done the my mind wondered and I started thinking about making pickups.

However,

John got the model working exactly the same as Trevors using a VBA app that avoids the rounding errors by storing results as txt files rather than numbers using the code from the book. From memory it was the last line in Dbar had the wrong sign compared to Trevors results but we could not find a mistake. The VBA coded seemed to solve the problem. Perhaps the rounding problem has the effect of changing signs when these small numbers are added and subtracted.

But happy there are more brains poking at this problem. I'll finish my code when my brain feels a bit clearer and it cools down a bit.

Cheers
Dom

[kiwigeo]

Take it to the propeller head's forum guys

[Jim Kirby]

But this IS the propeller head forum, isn't it? (Or at least the propeller head thread).

Dom, I got the same Dbar as the book, but different numerator in y_a. The step from AII 2-38 to AII 2-43 clearly has an error as printed, in the numerator.

Happy New Year everyone.

Jim

[Dominic]

Hey, you don't have to read if it makes your brain hurt Martin.

Jim, what do you think AII 2-38 to AII 2-43 should be?
I didn't mean the derivation of Dbar, I meant to result. Trevor posted his results for 100Hz (I think) and John and I were getting the wrong sign and/or number. But as I mentioned, John got the same result as Trevor with the right applitudes by not using excel or mathlab.

i wouldn't be so sure of myself to say the equations are clearly correct but I have had many eyes check them out as i am sure trevor did and we have not seen an error. And John's success confirms something. So I still think the error I am getting is the way excel calculates Dbar and so everything beyond that is going to be wrong (compared to Trevor's results).

Dom

[Steve]

Jim, what do you think AII 2-38 to AII 2-43 should be?

Hey Dom - I believe that it's the derivation of AII 2-38 to create AII 2-43. If you derive AII 2-38, we've both been getting F(+Kt.Kb.alphaab +...) instead of the F(-Kt.Kb.alphaab +...). A more complete explanation here: viewtopic.php?f=33&t=4636&view=unread&s ... 8e2#p57832.
Steve

[Jim Kirby]

Jim, what do you think AII 2-38 to AII 2-43 should be?

Hey Dom - I believe that it's the derivation of AII 2-38 to create AII 2-43. If you derive AII 2-38, we've both been getting F(+Kt.Kb.alphaab +...) instead of the F(-Kt.Kb.alphaab +...). A more complete explanation here: viewtopic.php?f=33&t=4636&view=unread&s ... 8e2#p57832.
Steve

Right, there should be a plus sign in front of the first term in the numerator in AII 2-43 instead of a minus sign. That follows directly from computing the determinant in AII 2-38.

[Dominic]

OK, I have not had the mental fortitude to go back and start trolling through all my equations again yet. But I did check the VBA chart that John created and I was wrong. The peaks are all correct but the amplitude of the air is a little low.
Again, from memory, when plugging Trevors 100Hz numbers for Dbar into the 4 Y equations I think we got the same answers as Trevor. But its too early in the year for my brain just yet. I'll follow on this thread and make any changes suggested but there is still an issue with the air res even when the sign is changed. So there is more work to do yet but the more eyes the better.

The important point is, Trevor got the right results, the model being a very good fit to the guitar. And that was using excel so I believe there must be a way to get around excel's shortcomings.

Cheers
Dom

[Jim Kirby]

Dom,

I'm curious about the mismatch in the air resonance peak after my change as well, given that the rest of the plot is really close to Trevor's. I want to spend some time checking on what parameters that peak is sensitive to - maybe something that isn't all that well constrained? (i.e., can I regain the fit to the experimental data by reasonable adjustments of material properties?) Is the spectral analysis of the experimental data subject to comparable computational uncertainties itself? (I think Trevor uses VA usually? I have no idea what the guts of VA's computational engine is based on.)

But right now, I'm too busy keeping certain parts of the world's infrastructure safe from tsunamis (or at least telling them in advance how deep a pile they are in if anything happens).

Jim

[Trevor Gore]

Thanks for keeping this going everyone, and sorry for not getting back. I took a holiday! Then there's all the other sh!t I've been having to deal with recently (nothing serious, just takes time). Anyway, there's about two days effort to go through all my workings again to check everything term by term. It's that mental fortitude thing... I'll get to it eventually...

Please keep posting.

BTW, has anyone tried yet to fit their version of the model to a guitar frequency response curve?

[Jim Kirby]

Thanks for keeping this going everyone, and sorry for not getting back. I took a holiday! Then there's all the other sh!t I've been having to deal with recently (nothing serious, just takes time). Anyway, there's about two days effort to go through all my workings again to check everything term by term. It's that mental fortitude thing... I'll get to it eventually...

Please keep posting.

BTW, has anyone tried yet to fit their version of the model to a guitar frequency response curve?

I've gotten too busy myself. In connection with your question, I'm wondering what parameter the main air peak is most sensitive to and whether that is possibly one of the more poorly constrained parameters. The "corrected" model still fits your plot pretty well except for the amplitude of the main air peak.

Direct answer to your question is no, not here. I don't have data that good for any of my builds so far, or any of them still sitting around to thump on.