Code rewrite!!

Hi again! I was searching the web for other people writing numerical solvers in Haskell and found this nice guy called James. His idea is to extract the tolerance test to a new function. This function takes an infinite list of approximated x values and then take the value that gives a relative error less than the specified tolerance.

Now relative error is calculated as absolute value of (x(n-1)-x(n))/x(n-1). I changed his convergence function and corrected a mistake in it, cause he forgot to take the absolute value for the whole expression of relative error. Here is my code:

module HaskNum where

converge :: (Ord a, Fractional a) => a -> [a] -> a
converge tol (old:new:xs) = if reldiff old new <= tol
                          then new
                          else converge tol (new:xs)
    where reldiff old new = abs $ (old-new) / old
converge tol (x:_) = x -- to handle when some functions
                       -- succeed from the first time

What I did was put it in a module I call HaskNum that will grow with time. And here is the new version of my previous four functions:

bisec :: (F -> F) -> (F,F) -> [F]
bisec f (a, b)
 | f(c) == 0 = [c]
 | f(a)*f(c) < 0 = c : bisec f (a,c)
 | otherwise = c : bisec f (c,b)
  where c = (a+b)/2 

fpi :: (F -> F) -> F -> [F]
fpi = iterate

newton :: (F -> F) -> (F -> F) -> F -> [F]
newton f f' = iterate xc
 where xc xi = xi - f(xi)/f'(xi)

secant :: (F -> F) -> (F, F) -> [F]
secant f (x0, x1) = xc : secant f (x1, xc)
 where xc = x1 - f(x1)*(x1-x0)/(f(x1)-f(x0))

Now things look really nice! And oh here how u can call one of them with a function that u define:

main = do
 let answer = (converge 0.0000005 . bisec (func)) (0,1)
 print answer

An interesting note is that all the methods except the Bisection Method gives the answer xc=0.6823278. The Bisection Method gives 0.682328, i.e. one digit less. Now that is a mystery for me, but I would guess its something with computer arithmetic and the intermediate operations.

Bisection Method vs Fixed-Point Iteration

How much faster is FPI compared to the Bisection Method? For that we need to think about accuracy. Or the error compared to the true solution. Now I know we don't know the true solution otherwise we would not be needing these numeric methods! But we can enclose the true solution inside an interval where we know that it lies in it (well for a continuous functions in the interval of interest).

Starting from The Bisection Method we know for each step the interval [a,b] is cut in half! After n steps our new interval [an, bn] have the length (b-a)/2^n. But choosing the solution to be the midpoint of the n-interval we know the true solution (if it is not the midpoint) is within an interval of length (b-a)/2^(n+1).

A solution is correct within p decimal places if the error |xc-r|<0.5*10^-p.

How many times do we need to cut the interval [0,1] in half to get an approximation with 6 correct digits for the problem f(x)=x^3+x-1 using the Bisection Method?

Solving (b-a)/2^(n+1) < 0.5*10^(-6) => n = 19.9. Which means The Bisection method will need 20 steps for 6 correct digits. The approximated solution with 6 correct digits is xc=0.682327.

Using FPI with g(x)=(1+2x^3)/(1+3x^2) and an initial guess x0=1. After 4 steps we get the same solution. A major improvement of number of steps by using FPI with a "nice" g(x)!

Fixed-Point Iteration

Fixed-Point Iteration, FPI for short, is another equation solving method. But! Unlike The Bisection Method you have to exercise more caution.

First we need to rewrite the initial problem from f(x)=0 to g(x)=x and iterate until the solution stabilizes. Here is the code for FPI:

fpi :: (F -> F) -> F -> Integer -> [F]
fpi _ _ 0 = []
fpi g x0 k = xc : fpi g xc (k-1)
 where xc = g(x0)

The code is simple and the method returns a list with the intermediate approximated solutions. The best solution is in the end of the list.

Let us try to solve the same equation as the one presented in The Bisection Method, f(x)=x^3+x-1. Re-writing it we would get the equation: x=1-x^3 thus g(x)=1-x^3 (remember it should be of the form g(x)=x). But FPI will fail and the list returned would look like:

[0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,...]

That is no good! But if we happen to add 2x^3 to both sides of the equation, and rearrange things we get: x = (1+2x^3)/(1+3x^2). So using g(x)=(1+2x^3)/(1+3x^2) we will have a good approximation in just 6 steps!

[0.75,0.68604654,0.68233955,0.68232775,0.6823278,
0.6823278]

Why? Well because there is a theorem that says: Assume g(x) is continuously differentiable, and g(r)=r, and |g'(r)|<1. Then FPI converges to the fixed point r for sufficiently close initial guess to r.

Checking |g'(r)| (r≃0.68) for the first attempt we can see that it is larger than 1, but not the second attempt.

FPI can converge to a solution faster than The Bisection Method, but we need to re-write the original solution. Because we dont now the solution r beforehand we can test re-writing it differently until FPI succeeds.