[SOLVED] Solving trigonometric equation

jfortiv · #1 $Old$ November 5th, 2009, 08:47 PM

Hi,

I'm trying to solve a simple equation in MATLAB, numerically:

tan(x) = tanh(x)

I'm not very skilled in MATLAB and I've poked around in the docs for 15 minutes or so now. I looked at solve() in the symbolic toolbox. I can't get that to give me anything other than zero. I basically want the first 20 or so non-zero solutions for x (there is an infinity of solutions to this equation).

Any guidance (even just pointing me to the right docs) would be greatly appreciated.

Thanks!

elbarto · #2 $Old$ November 5th, 2009, 10:52 PM

This will find a root to the equation numerically by using an inital trial solution of x = pi:

Code:

EDU>> S = @(x) tanh(x)-tan(x)

S = 

    @(x)tanh(x)-tan(x)

EDU>> fsolve(S,pi)
Optimization terminated: first-order optimality is less than options.TolFun.

ans =

    3.9266

EDU>>

Not sure on how to get 20 solutions numerically unless you can find values close the the actual solutions for each case.

Elbarto

CaptainBlack · #3 $Old$ November 6th, 2009, 03:40 AM

Quote:

Originally Posted by jfortiv $View Post$

Hi,

I'm trying to solve a simple equation in MATLAB, numerically:

tan(x) = tanh(x)

I'm not very skilled in MATLAB and I've poked around in the docs for 15 minutes or so now. I looked at solve() in the symbolic toolbox. I can't get that to give me anything other than zero. I basically want the first 20 or so non-zero solutions for x (there is an infinity of solutions to this equation).

Any guidance (even just pointing me to the right docs) would be greatly appreciated.

Thanks!

Rewrite the equation in the form: $\tan(x)-\tanh(x)=0$ , then there is one root of this in each of the intervals $(\pi(n-1),\pi(n+1)),\ \ \ n=0,\pm 1,\pm 2, ...$ .

Now start your solver at the points $2n\pi$ for a suitable selection of $n$ 's to get your 20 solutions.

CB

jfortiv · #4 $Old$ November 6th, 2009, 05:42 AM

Elbarto and CaptainBlack,

Thank you so much.

Can you show how you obtained these intervals?

I'm just going to solve using a for loop with an input vector of xo values. It's seems a shame that this is in a special MATLAB toolbox. I would think this kind of solve would be core functionality.

elbarto · #5 $Old$ November 6th, 2009, 07:13 AM

jfortiv,

You can use "fzero" instead of "fsolve". I use "fsolve" as a force of habit. "fzero" will be fine for this case I believe. The difference between the two is that "fzero" will only find a root if the function cuts the x-axis, whereas "fzero" will find the root if the function touches the axis or cuts it. I may be wrong tho so someone please correct me if this is incorrect.

As for the intervals, I think this is ok (although the initial guess chosen in the range that CB described appears to be the actual roots?)

Code:

n = 1:20;%number of roots
B = [pi.*(2.*n-1)' pi.*(2.*n+1)'];%bounds CB described
f_roots = zeros(length(n),1);%preallocate memory
f = @(x) tanh(x)-tan(x);%nonlinear equation to solve
for i = n
    start_point = mean(B(i,:));%mid point of bounds
    f_roots(i) = fzero(f,start_point);%solve root
end

Elbarto

CaptainBlack · #6 $Old$ November 6th, 2009, 11:15 AM

Quote:

Originally Posted by jfortiv $View Post$

Elbarto and CaptainBlack,

Thank you so much.

Can you show how you obtained these intervals?

$\tan(x)$ is periodic with period $\pi$ , is continuous in each of the given open intervals and goes from $-\infty$ to $\infty$ in each of these intervals, so there is at least one root in each interval. Now look at the form of $\tanh(x)$ and it is obvious that there is but one root in each of the intervals.

CB

jfortiv · #7 $Old$ November 6th, 2009, 11:43 PM

This is interesting...

I did a bit of experimenting and found that if I start with an initial guess near the asymptote of tan(x), fsolve() will pick up the asymptote as a zero. So I need to start with initial guesses near the zeros that I want.

Maybe I'm being naive, but I find this frustrating because it requires some intuition about the solution before obtaining the solution values numerically. I would think that one should be able to tell MATLAB, "solve tanh(x)=tan(x) for all zeros in the range 0<x<=20, and don't give me any asymptotes as a solution." Is this technology just not available?

Note: I needed to divide the end points of the intervals by 1/2 because I noticed that I was skipping roots. I think this is right because the solutions are spaced out by about pi if you plot tanh(x) - tan(x). Am I right in concluding that the choice of intervals is merely an observational process?

Thanks Elbarto and CaptainBlack. I'm learning a lot by working through this problem.

CaptainBlack · #8 $Old$ November 7th, 2009, 12:04 AM

Quote:

Originally Posted by jfortiv $View Post$

Note: I needed to divide the end points of the intervals by 1/2 because I noticed that I was skipping roots. I think this is right because the solutions are spaced out by about pi if you plot tanh(x) - tan(x). Am I right in concluding that the choice of intervals is merely an observational process?

That's a mistake on my part, $\tan(x)$ is periodic with period $\pi$ not $2\pi$ (my earlier posts have been corrected). Also when you plot the functions it is obvious that the roots all lie near the mid point of the intervals between the infinities of $\tan$ .

CB

CaptainBlack · #9 $Old$ November 7th, 2009, 12:47 AM

Quote:

Originally Posted by jfortiv $View Post$

This is interesting...

I did a bit of experimenting and found that if I start with an initial guess near the asymptote of tan(x), fsolve() will pick up the asymptote as a zero. So I need to start with initial guesses near the zeros that I want.

Maybe I'm being naive, but I find this frustrating because it requires some intuition about the solution before obtaining the solution values numerically. I would think that one should be able to tell MATLAB, "solve tanh(x)=tan(x) for all zeros in the range 0<x<=20, and don't give me any asymptotes as a solution." Is this technology just not available?

This is one of the things I don't like about Matlab, it appears to give the impression that numerical analysis is a solved problem and that there can be a tool for each numerical question that can be asked.

For one example near the singularities of $\tan$ the underlying assumptions of fsolve do not hold.

For a realistic view of numerical analysis you should read Forman Acton's "Numerical Methods that usually work"

Amazon.com: Numerical Methods that Work (Spectrum) (9780883854501): Forman S. Acton: Books

(if you read the documentation on fsolve and Acton's book something will strike you as odd. The basic idea in fsolve appears to be finding a root of $f(x)=0$ by minimising $(f(x))^2$ , this is a pet hate of Acton's and he will explain why, in general, this is a bad idea)

CB

elbarto · #10 $Old$ November 7th, 2009, 01:49 AM

Quote:

Originally Posted by jfortiv $View Post$

This is interesting...

Maybe I'm being naive, but I find this frustrating because it requires some intuition about the solution before obtaining the solution values numerically. I would think that one should be able to tell MATLAB, "solve tanh(x)=tan(x) for all zeros in the range 0<x<=20, and don't give me any asymptotes as a solution." Is this technology just not available?

I can understand your frustration as there are somethings that MATLAB doesn't provide out of the box that seem like they should be there. As far as I am aware there is no function that will allow you to find multiple roots of an arbitrary periodic function.

That said, if you understand the basics of a few different methods of root finding, you might be able to use something like bisection method to solve the problem without knowing to much about the function behavior (although this is a rather ugly and inefficient approach I think)

Here is some code for the Bisection or Bracketing method:

Code:

function [x fval iter] = mhfBracket(f,xi,xj,MAXERR,MAXIT)
% f      -   function_handle
% xi     - left bound - double
% xj     - left bound - double
% MAXERR - Maximum Error - double
% MAXIT  - Max Iterations - int
if nargin < 4;MAXERR = 1e-6;end
if nargin < 5;MAXIT = 100;end
iter = 0; %start counter
while abs(xj-xi) > 2*MAXERR && iter < MAXIT ;
    iter = iter +1; %keep track of iterations    
    xm = (xi+xj)/2; %mid point
    if f(xi)*f(xm) > 0    
        xi = xm;
    else
        xj = xm;
    end
end
x = (xi+xj)/2;
fval = f(x);

Using this code to find the roots, and iterative process can be establish to find the bounds which bracket a root. In the code below I specify a value called "step" which I use to progressively increase my guess for an upper bound until the solution in bracketed.

Code:

function [froots fval iter] = mhfPeriodicBracket(f,xi,nroots,step,MAXERR,MAXIT)
if nargin < 5;MAXERR = 1e-6;end
if nargin < 6;MAXIT = 100;end
froots = zeros(nroots,1);
fval = zeros(nroots,1);
iter = zeros(nroots,1);
for r = 1:nroots
    xj = xi + step;%first xj guess
    while f(xi)*f(xj) > 0%look for suitable xj
        xj = xj + step;
    end
    [x e i] = mhfBracket(f,xi,xj,MAXERR,MAXIT);
    froots(r) = x;
    fval(r) = f(x);
    iter(r) = i;
    xi = x + 2*MAXERR;% get new xi value
end

Using this function I can now find the first 20 roots starting with with an initial guess at 0:

Code:

EDU>> [froots fval num_iterations] = mhfPeriodicBracket(f,0,20,pi/10)

froots =

    0.0000
    1.5708
    3.9266
    4.7124
    7.0686
    7.8540
   10.2102
   10.9956
   13.3518
   14.1372
   16.4934
   17.2788
   19.6350
   20.4204
   22.7765
   23.5619
   25.9181
   26.7035
   29.0597
   29.8451


fval =

  1.0e+008 *

   -0.0000
    0.0284
    0.0000
   -0.0135
   -0.0000
   -0.0157
    0.0000
   -0.0188
    0.0000
   -0.0235
    0.0000
   -0.0312
    0.0000
   -0.0464
    0.0000
   -0.0906
   -0.0000
   -1.9091
   -0.0000
    0.1001


num_iterations =

    18
    20
    21
    19
    21
    19
    21
    19
    21
    19
    21
    19
    21
    19
    21
    19
    21
    19
    21
    19

EDU>>

This method may or may not work for what you want to do, you can always modify the criteria to quit if a root is greater then a certain value rather then specifying how many roots to find like I did. The problem with this approach is that I am also picking up the asymptotes as roots but you can modify the code accordingly (ie continue next search if fval > some_value). This is a very basic numerical approach and while it might work ok with simple functions like sin(x), since your function has asymptotes you need to be careful in trusting the output but it seems like it might be useful when only little information is known about the function.

Just my 2c attempt anyway mate -

Regards Elbarto

elbarto · #11 $Old$ November 7th, 2009, 01:54 AM

CB,

What level of reading is this book at? ie is the content able to be comprehended with only a general engineering background? I wouldn't mind learning a little more about numerical methods since I take them for granted quite a bit.

Regards Elbarto

CaptainBlack · #12 $Old$ November 7th, 2009, 03:23 AM

Quote:

Originally Posted by elbarto $View Post$

CB,

What level of reading is this book at? ie is the content able to be comprehended with only a general engineering background? I wouldn't mind learning a little more about numerical methods since I take them for granted quite a bit.

Regards Elbarto

If your mathematical background is such that you have no difficulty with Kreysig's Engineering Maths, then that is more than enough to read Acton.

CB

CaptainBlack · #13 $Old$ November 7th, 2009, 03:26 AM

Quote:

Originally Posted by elbarto $View Post$

That said, if you understand the basics of a few different methods of root finding, you might be able to use something like bisection method to solve the problem without knowing to much about the function behavior (although this is a rather ugly and inefficient approach I think)

Bisection is usually my first choice, but you have to remember that it is not a root finder, but a zero crossing finder, which is often adequate for finding a root, but it will not find a root of even multiplicity.

CB

jfortiv · #14 $Old$ November 7th, 2009, 02:09 PM

I'm not familiar with the usage of these methods, but I'm very interested. Thanks for providing the code and the reference. I'll research this stuff a bit when I find some time.

My friend also has a book called "Analysis of Numerical Methods" by Eugene Isaacson and Herbert Bishop Keller. I've only briefly glanced at their description of the Newton-Raphson method.

If either of you have any interest, I asked about solving this equation because I wanted to solve for the eigenvalues of the spatial portion of the wave equation for an Euler-Bernoulli beam clamped at one end and pinned at the other. I put together an M-file (attached as dynBeam.txt) that solves for a steel beam over time and animates the motion with some splined curve that I threw together for the IC.

Thanks again!