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?
Walter Roberson
Walter Roberson 2017년 4월 24일
It is not clear what your question is?

댓글을 달려면 로그인하십시오.

답변 (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 CenterFile Exchange에서 Numeric Solvers에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by