Fixed point iteration to find the roots of 0=x-tan(x)

[Brendy]

Homework Statement

The question wants me to first estimate the roots by drawing the graph and then by using a 'suitable' fixed point method to determine the first 4 positive roots.

Homework Equations

0=x-tan (x)
I rearranged to get x=arctan (x) so that the series x_n will converge.

The Attempt at a Solution

I've attached my code. My lecturer told us to reuse some code in the lecture notes and just modify it for this function but I don't fully understand what it's doing. My main problem is the f_old and f_new parts. That was my attempt at changing what was in the book (the book had a different function so I had to modify this for my function) but I don't know how they're useful.
Anyway, if I let the code run, my x_new and f_new go to zero which I think means that it's honing in on the x=0 root. The next root is about x=4.5 or so but if you add pi to 0 then you get pi and not the next root (x=4.5). What am I missing?

 

[vela]

Did you graph tan x and x and see where the two intersect? It's pretty clear from that what's going on.

 

[Brendy]

I did graph them. From estimating by eye, the roots were: x=0, 4.5, 7.8 and 11 +/- 0.5. What's clear? I'm pretty new to MATLAB so excuse my incompetence in it.

 

[vela]

Why would you expect the next root to be at x=π then?

 

[Brendy]

Wouldn't adding $$\pi$$ to the root give you the next root because tan (x) repeats?

 

[vela]

No, the function y=x-tan x isn't periodic like tan x alone is, so the roots don't occur at constant intervals.

Did you graph the function y=x-tan x alone, or did you graph the two functions y=x and y=tan x and see where they intersect? If you do the latter, I think you'll see why the roots get farther and farther apart as x increases.

 

[Brendy]

We were told to graph y=x and y=tan x on the same graph and find where they intersect. Ah, the roots are all along y=x so adding pi would only get the roots if it was y=c, correct? How would I get the next root though? I thought maybe adding root(2) *pi would get me the next root but it doesn't.

 

[vela]

Yes, the roots get farther apart because y=x increases as you move right to the next root. There's no simple pattern to find the exact roots, but as x gets large, the line y=x hits tan x where it also gets large. Where would you say that approximately happens?

 

[Brendy]

I'm sorry but I don't understand what you're asking. Do you mean for what values of y the y=x line intersects tan x? If so, then the y=x line hits tan x whenever x=tan x. The x and y values are the same. Is that what you were asking?

 

[vela]

No, I was just saying that as x gets large, tan x must also be large (because you want x=tan x). This means the roots are near the odd multiples of pi/2, where tan x blows up.

 

[Brendy]

Oh right, close to the asymptotes. That doesn't really help actually finding the exact roots though. My lecturer must have been mistaken when he told us how to find them. He was saying that when we run our program, we'll find a root between -pi/2 and pi/2 and to get the next one we'd have to go atan (x) + n*pi to get the next one. I don't get the roots if I do that though.

 

[vela]

He might have just meant that you can shift by pi to come up with an initial guess for the next root. Usually, when you find roots numerically, you need to give the algorithm a starting point. Depending on which point you start at, the algorithm will converge on different roots.

 

[Brendy]

Of course! Thanks, that helps a lot. How would I put that into the code though? I was thinking it would be a for loop of some sort outside of the loop that's doing the rest of the calculations with n=0:4 but I'm not sure how it would look. Here's my attmept:

for n=0:4
---while looping %start iterations
------x_new=atan(x_old);
------f_old=atan(x_old)-x_old;
------f_new=atan(x_new)-x_new;

------loop=loop+1;
------if loop>loop_max
---------looping=false;
------end

------if abs(x_new-x_old)<small_number
---------looping=false;
------end

------if abs(f_new)<small_number
---------looping=false;
------end

------x_old=x_new;
------disp([loop x_new f_new]);
---end
x_old=atan(x_old)+n*pi;
end

I tried it and it only gives me one table of values which are converging to x=0.

 

[vela]

You don't want to add n*pi; you want to just add pi to the previous root. You also need to reset looping to true before reentering the while loop.

 

[Brendy]

How would I do that? Would just adding the line, looping=true at the end of the while loop do it?

 

[vela]

Sounds reasonable. Try it out.

 

[Brendy]

I added just the looping=true at the end of the while loop and that only made it do another 4 iterations so I added looping=0 so it would repeat the while loop and do 50 iterations again. For some reason it didn't do it. It did 55 iterations and stopped. Here's what my code looks like now:

x_old=1;----------%initial guess
loop=0;
loop_max=50;----------%max number of iterations
small_number=0.00001;----------%target accuracy
looping=(loop<loop_max);----------%looping as long as loop<loop_max

for n=0:4
------while looping----------%start iterations
------x_new=atan(x_old);
------f_old=atan(x_old)-x_old;
------f_new=atan(x_new)-x_new;

------loop=loop+1;
------if loop>loop_max
----------looping=false;
------end

------if abs(x_new-x_old)<small_number
----------looping=false;
------end

------if abs(f_new)<small_number
----------looping=false;
------end

------x_old=x_new;
------disp([loop x_new f_new]);
---end
---looping=0;
---looping=true;
---x_old=atan(x_old)+pi;
end

 

[Brendy]

Also, it seems that no matter what I put as the initial guess, the algorithm converges to the x=0 root. Is there something I've forgotten to add to my code?

 

[D H]

Of course that is only going to find the root at x=0. That code is solving atan(x)=x, which is *not* the same as x=tan(x). The former problem has one solution; the latter has an infinite number of solutions.

 

[Brendy]

You're right. My lecturer was talking about using atan (x)=x instead of tan(x)=x because the derivative of atan(x) converges while the derivative of tan(x) does not. Can you explain that relative to the code? I'm incredibly confused and the lecture notes and textbook aren't helping.

 

[D H]

The derivative of tan(x) does not diverge. It is quite well defined. The problem is that tan'(x)=sec²x >= 1 for all x. This guarantees that a fixed point iteration with x_n+1=tan(x_n) will not (cannot!) converge. Suppose we are trying to find the solution just below 10.5*pi. The true solution is 32.95638904... Starting with an initial value of 32.95640764 (a very good guess, and easy to generate), the sequence generated from x_n+1=tan(x_n) is (32.95640764, 32.97661717, 98.9507368, 106.0213008, -1.015014699, ...) So even though the initial guess was very, very good, this method fails.

The derivative of atan(x) is 1/(1+x²), so 0 < atan'(x) <= 1. Except for those points where tan(x)=0, that <= becomes < -- and this is the true condition that guarantees that a fixed point iteration will converge. BTW, that atan'(x)=1 at x=0 means that using a fixed point iteration to find the solution at x=0 will converge very, very slowly. All of the other solutions will converge quite rapidly if you can come up with an atan-based solution.

Of course to do that you need to get around the fact than -pi/2 < atan(x) < pi/2. You will need to rewrite the problem a bit using trig identities to accomplish that.

 

[Brendy]

Ok, I think I follow. How would I rewrite the problem? Do you mean by substituting x=atan(x) into 0=x-tan (x)? That would give me atan (x) - x=0 which is atan(x)=x. It doesn't really get me anywhere.

 

[D H]

Try finding a new variable that pertains just to one particular solution. You know that the solutions for positive x fall between n*pi and n*pi+pi/2 for positive n. Take advantage of this.

 

[Brendy]

Ok, I don't follow anymore. tan (x) will give me a whole lot of solutions, each separated by pi, atan (x) only has one solution... I'm not sure what to do from there. How would I take advantage of x falling between n*pi and n*pi/2?

 

[D H]

Ok, I don't follow anymore. tan (x) will give me a whole lot of solutions, each separated by pi, atan (x) only has one solution... I'm not sure what to do from there. How would I take advantage of x falling between n*pi and n*pi/2?

First off, the solutions to x=tan(x) are not separated by pi. The separation asymptotically approaches pi. How do you take advantage of x falling between n*pi and (n+1/2)*pi? Write x as sum of two entities. Learn to invent new variables.

 

[Brendy]

So x=n*pi + theta. How do I determine what theta is? Theta<pi/2 but only by a little and it approaches pi/2 as n approaches infinity. I don't know how to find an exact solution though.

 

[D H]

You aren't going to find an exact solution. You need to use some estimation algorithm -- such as a fixed point iteration algorithm. With x=n*pi + theta, what does x=tan(x) say about theta?

 

[Brendy]

Since tan(n*pi)=0 it says that x=tan (theta)? Theta=arctan (x). So tan(theta)=n*pi +arctan(x). I'm I veering off course? It kind of feels like it :P

 

[D H]

Since tan(n*pi)=0 it says that x=tan (theta)?

Correct.

Theta=arctan (x). So tan(theta)=n*pi +arctan(x). I'm I veering off course? It kind of feels like it :P

Veering off course. Make the substitution on both sides of the equality, and remember that a fixed point iteration based on tangent is not stable. You need to use atan.

Don't worry about x until you have found theta.

 

[Brendy]

so n*pi<tan(theta)<(n+1/2)pi. So, atan(n*pi) < theta < atan((n+1/2)pi) and atan(n*pi)=~1.263, atan((n+1/2)pi=~1.362 so 1.263<theta<1.362.
That narrows it down a fair bit. I'm not sure what the exact forms of atan(n*pi) and atan((n+1/2)pi) are though.