Can't find all roots to function with fixed point method

Question

I have a function $f(x) = x^2-8x-10\cos(2x)+15$ and I'm supposed to via Matlab find all roots for this function. I can see graphically that it has $4$ roots and I've written a matlab function to find a root. When the absolute difference between $x_n$ and $x_{n+1}$ is smaller than $10^{-10}$ the function returns $x_{n+1}$ as the root.

Since it has $4$ roots I call the function $4$ times, with different starting values each time, all of them being close to some root.

The problem is this: the 1st and 3rd root are found. And the 2nd and 4th roots aren't. Trying to find the 2nd root returns me the value of the 1st root, and trying to find the 4th root just loops telling me $x_n$ is NaN.

Even when the start value when I try to find the 2nd root is the root, it still returns the 1st root.

This is for a school lab and it has to be found with a fix-point Method I've been given:

x_n = (1/15)x^2+(7/15)x-(2/3)cos2x+1

My function for finding roots looks like this

function root = fixpoint(x, T)

while true
    x_n = 1/15*(x^2)+7/15*x-2/3*cos(2*x)+1;
    error = abs(x_n-x);
    disp(['Error: ', num2str(error),' X_n: ', num2str(x_n)])

    if error < T
        break
    end

    x = x_n;

end

root = x_n;

end

And I call it 4 times with:

errorMargin = 10^-10;
p1 = fixpoint(2, errorMargin);
p2 = fixpoint(3.977, errorMargin);
p3 = fixpoint(5, errorMargin);
p4 = fixpoint(7, errorMargin);  

And I get: p1= 2.429(correct), p2= 2.429(actual is 3.977), p3= 5.571(correct), p4 loops 'NaN'

What do I need to do to find all $4$ roots?

Answer 1

Your problem is that fixed point iteration only converges when the derivative of the right hand side is less than $1$ in absolute value. If you linearize the problem, the error is multiplied by the derivative of the RHS at each step. If that derivative is greater than $1$ in absolute value, you will get farther and farther from the root. The value of the derivative at $3.977$ is about $2.3$ You don't say what the fourth root is, but you started at $7$ and the derivative looks to be greater than $1$ there as well. A plot of the derivative from Alpha is below.

Answer 2

You ar looking for fixpoints of $g(x)=c\cdot f(x)-x$. It is not surprising that two of the four fixpoints intersect the diagonal $y=x$ from below, i.e., at a slope $g'(x)>1$. Such fixpoints are repelling and will take you further away with each iteration. Only fixpoints where $|g'(x)|<1$ can be found by iterating from a sufficientla close starting point.

The solution is to try different $c$, not only (as you did) $c=\frac1{15}$. What change to $c$ could make the graph of $g$ make the diagonal intersect from above instead?

Answer 3

That's the problem with fixed point iteration: $$g(r+\epsilon)\approx g(r)+\epsilon g^{\prime}(r)$$ So that if $\left|g^{\prime}(r)\right|>1$, successive iterations will grow away from the root. You were prescribed a fixed point iteration to use, so how do you make it work in this face of this adversity? A little analysis of the problem can help: if $x_0=r+\epsilon$, then $x_1\approx r+\epsilon g^{\prime}(r)$ and $x_2\approx r+\epsilon \left(g^{\prime}(r)\right)^2$. Then $$\frac{\left|x_2-x_1\right|}{\left|x_1-x_0\right|}\approx\left|g^{\prime}(r)\right|$$ And we know that if $\left|g^{\prime}(r)\right|>1$ we are in trouble. But in this case, $$\epsilon\approx\frac{\left(x_1-x_0\right)^2}{x_2-2x_1+x_0}$$ So we can subtract this from $x_0$ and hope to approach $r$. Here is my matlab function:

% fixpoint.m

function x2 = fixpoint(f,x0);

quad = 0;
linear = 0;
x1 = f(x0);
err = 1;
calls = 1;
while abs(err) > 1.0e-10,
    x2 = f(x1);
    calls = calls+1;
    if abs(x2-x1) > abs(x1-x0),
        err = (x1-x0)^2/(x2-2*x1+x0);
        x0 = x0-err;
        x1 = f(x0);
        calls = calls+1;
        quad = 1;
    else
        err = x2-x1;
        x0 = x1;
        x1 = x2;
        linear = 1;
    end
end
x2 = x0;
if quad,
    fprintf('Quadratic method was used\n');
end
if linear,
    fprintf('Linear method was used\n');
end
fprintf('calls = %d\n',calls);

And my script that drives it:

% fixtest.m

f = @(x) x.^2/15+7/15*x-2/3*cos(2*x)+1;
for x = [2.4 3.9 5.4 6.8],
    y = fixpoint(f,x);
    fprintf('%f %13.11f\n',x,y);
end

And the output:

>> fixtest
Linear method was used
calls = 33
2.400000 2.42974456196
Quadratic method was used
calls = 9
3.900000 3.97704805536
Linear method was used
calls = 11
5.400000 5.57154144613
Quadratic method was used
Linear method was used
calls = 10
6.800000 6.72106227846

Note that fixpoint.m always detected whether $\left|g^{\prime}(r)\right|>1$ correctly right off the bat so there was never any switching between the linear and quadratic methods. Now, if $g^{\prime}(r)=1$ I would still have been toast...