Nonlinear Tangent (trigonemetric) equation

Could you guys help me the solve following problem. This is the dispersion equation for a symmetric lamb wave. I am looking for the z values for different f (frequencies)? I have attached a paper which has the original equation (equation #1 ). I have to find something similar to Fig 1 a as in the attached folder.

 close all

 clear all

 clc

 %cl longitudunal wave speed

 %ct transversal wave speed

 %cp phase velocity of wave

 %cg group velocity of wave

 % z=ct/cp;

 pi=3.14;

 nu=0.33;% poisson ratio

 ro=2700;%density kg/m3

 E=70e9;% elastic modulus Pa

 mu=E/(2*(1+nu)); %shear modulus

 cl=((E*(1-nu))/(ro*(1+nu)*(1-2*nu)))^(0.5);

 ct=(mu/ro).^(0.5);

 k=ct/cl;

 f=10:10:3e6;

 w=2*pi*f;

 d=(w*0.8e-3)/ct;

 fzero (@(z) (2*z.^2-1).^2*(sin(sqrt((1-z.^2))*d))*cos(sqrt((k.^2-z.^2.*d)))-(sin(sqrt  (k.^2-z.^2.*d)))*cos(sqrt((1-z.^2))*d)*(4*z.^2)*sqrt(1-z.^2)*sqrt(k.^2-z.^2),1)

[Andrew Newell]

It's really hard to see how your code relates to Equation 1 in the test. Rather than forcing us to guess, could you harmonize the terminology? For example, is z the same as h in Equation 1? And where did alpha, beta and the tanh functions go?

[Walter Roberson]

It is not clear what your question is?

Answer 1

[Roger Stafford]

You are not using 'fzero' correctly. It cannot use a function handle that produces a vector as yours does - in addition to a changing 'z', for different values of 'f' you have different values of 'd'. You will have to call on 'fzero' for each different value of 'd' in a for-loop, and that will be 300,000 calls! It might pay you to find a more efficient starting value than 1 so as to speed up each call, for which some plots of your function could be of considerable help in predicting the solutions.

Answer 2

[cagatay yilmaz]

Dear all

My apologize to not provide enough argument. I go through the final equation and revised it as indicated below. You can find the my steps that shows how I get to final equation in the attached pictures. The final equation can be found also last line of 2.jpg. I explained what happened to alpha, beta, k and other variables in the attached pictures. Now when i run the script it works but it produces imaginary results.

 close all

 clear all

 clc

 %cl longitudunal wave speed

 %ct transversal wave speed

 %cp phase velocity of wave

 %cg group velocity of wave

 pi=3.14;

 nu=0.33;% poisson ratio

 ro=2700;%density kg/m3

 E=70e9;% elastic modulus Pa

 mu=E/(2*(1+nu)); %shear modulus

 cl=((E*(1-nu))/(ro*(1+nu)*(1-2*nu)))^(0.5);

 ct=(mu/ro).^(0.5);

 h=1e-3; % half thickness

 for f=100:100:2e6
%      in(f/100)=5000-f/100;

 dt(f/100)=2*pi*f*h/ct;

 dl(f/100)=2*pi*f*h/cl;

sol(f/100)=fzero(@(cp) tan(dt(f/100)*sqrt(1-ct^2/cp^2))/tan(dl(f/100)*sqrt(1-cl^2/cp^2))+(4*sqrt(1-cl^2/cp^2)*sqrt(1-ct^2/cp^2)*ct^3)/((2*ct^2/cp^2-1)^2*cl*cp^2),-1000);
end

Answer 3

[Andrew Newell]

Before trying to find zeros for a function, it's a good idea to plot it so you understand the nature of the problem. If you do that with this function for f=100, you'll find that it is imaginary except for a particular range of cp and that it is always nonnegative. It touches the cp axis at a couple of points. You'll never find those zeros with fzero.

The good news is, the answer is actually very simple. If you have the Symbolic Toolbox, try the following:

syms dt ct cp dl cl
fun = tan(dt*sqrt(1-ct^2./cp.^2))./tan(dl*sqrt(1-cl^2./cp.^2))+(4*sqrt(1-ct^2./cp.^2).^1.5*ct^3)./((2*ct^2./cp.^2-1).^2*cl.*cp.^2);
pretty(fun)

    /        /       2 \ \        /       2 \3/4
   |        |     ct  | |      3 |     ct  |
tan| dt sqrt| 1 - --- | |    ct  | 1 - --- |    4
   |        |       2 | |        |       2 |
   \        \     cp  / /        \     cp  /
------------------------- + ---------------------
   /        /       2 \ \          /     2     \2
   |        |     cl  | |        2 | 2 ct      |
tan| dl sqrt| 1 - --- | |   cl cp  | ----- - 1 |
   |        |       2 | |          |    2      |
   \        \     cp  / /          \  cp       /

It should be obvious where the zeros are, and you can confirm it with:

subs(fun,cp,ct)

ans =

0

subs(fun,cp,-ct)

ans =

0

Note also that it goes to infinity at cp = +-cl and at cp = +- ct*sqrt(2).