Nonlinear Tangent (trigonemetric) equation

2 次查看（过去 30 天）
cagatay yilmaz
cagatay yilmaz 2017-4-23
回答： Andrew Newell 2017-4-25
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 个评论
Andrew Newell
Andrew Newell 2017-4-24
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

Roger Stafford
Roger Stafford 2017-4-24
编辑：Roger Stafford 2017-4-24
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 2017-4-25
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 2017-4-25
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).

类别

在 Help Center 和 File Exchange 中查找有关 Numeric Solvers 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by