Nonlinear Tangent (trigonemetric) equation

2 views (last 30 days)
cagatay yilmaz
cagatay yilmaz on 23 Apr 2017
Answered: Andrew Newell on 25 Apr 2017
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)
  2 Comments
Walter Roberson
Walter Roberson on 24 Apr 2017
It is not clear what your question is?

Sign in to comment.

Answers (3)

Roger Stafford
Roger Stafford on 24 Apr 2017
Edited: Roger Stafford on 24 Apr 2017
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.

cagatay yilmaz
cagatay yilmaz on 25 Apr 2017
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

Andrew Newell
Andrew Newell on 25 Apr 2017
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).

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!