# Fixed point iteration¶

Let’s convert the roots of a quadratic polynomial \(f(x)\) to a fixed point problem.

```
using FundamentalsNumericalComputation
```

```
p = Polynomial([3.5,-4,1])
r = roots(p)
@show rmin,rmax = sort(r);
```

```
(rmin, rmax) = sort(r) = [1.2928932188134525, 2.7071067811865475]
```

We define \(g(x)=x-p(x)\). Intersections of its graph with the line \(y=x\) are fixed points of \(g\) and thus roots of \(p\). (Only one is shown in the chosen plot range.)

```
g = x -> x - p(x)
plt = plot([g x->x],2,3,l=2,label=["\$y=g(x)\$" "\$y=x\$"],aspect_ratio=1,
title="Finding a fixed point",legend=:bottomright)
```

If we evalaute \(g(2.1)\), we get a value of almost 2.6.

```
x = 2.1; y = g(x)
```

```
2.59
```

So \(g(x)\) is considerably closer to a fixed point than \(x\) was. The value \(y=g(x)\) ought to become our new \(x\) value! Changing the \(x\) coordinate in this way is the same as following a horizontal line over to the graph of \(y=x\).

```
plot!([x,y],[y,y],label="",arrow=true)
```

Now we can compute a new value for \(y\). We leave \(x\) alone here, so we travel along a vertical line to the graph of \(g\).

```
x = y; y = g(x)
plot!([x,x],[x,y],label="",arrow=true)
```

You see that we are in a position to repeat these steps as often as we like. Let’s apply them a few times and see the result.

```
for k = 1:5
plot!([x,y],[y,y],label=""); x = y; # y --> new x
y = g(x); plot!([x,x],[x,y],label="") # g(x) --> new y
end
display(plt)
```

The process spirals in beautifully toward the fixed point we seek. Our last estimate has almost 4 accurate digits.

```
abs(y-rmax)/rmax
```

```
0.0001653094344995643
```

Now let’s try to find the other fixed point \(\approx 1.29\) in the same way. We’ll use 1.3 as a starting approximation.

```
plt = plot([g x->x],1,2,l=2,label=["\$y=g(x)\$" "\$y=x\$"],aspect_ratio=1,
title="No convergence",legend=:bottomright)
x = 1.3; y = g(x);
for k = 1:5
plot!([x,y],[y,y],label="",arrow=true); x = y; # y --> new x
y = g(x); plot!([x,x],[x,y],label="",arrow=true) # g(x) --> new y
end
display(plt)
```

This time, the iteration is pushing us *away* from the correct answer.