Skip to main content

Nicolás' blog

Numerical methods challenge: Day 2

Nicolás Guarín-Zapata

Comments

Also available in:

Español

During October (2017) I will write a program per day for some well-known numerical methods in both Python and Julia. It is intended to be an exercise then don't expect the code to be good enough for real use. Also, I should mention that I have almost no experience with Julia, so it probably won't be idiomatic Julia but more Python-like Julia.

Regula falsi

The second method to be considered is the false position method, or regula falsi. This method is used to solve the equation \(f(x) = 0\) for \(x\) real, and \(f\) continuous. It starts with an interval \([a,b]\), where \(f(a)\) and \(f(b)\) should have opposite signs. The method is similar to bisection method but instead of halving the original interval it takes as a new point of the interval the intercept of the line that connect the function evaluated at the two extreme points. Then, the new point is computed from

\begin{equation*} c = \frac{a f(b) - b f(a)}{f(b) - f(a)} \end{equation*}

As in bisection method, we keep the interval where the solution appears, i.e., the sign of the function changes.

We will use the function \(f(x) = \cos(x) - x\) to test the codes, and the initial interval is \([0.5, \pi/4]\).

Following are the codes.

Python

from __future__ import division, print_function
from numpy import abs, cos, pi

def regula_falsi(fun, a, b, niter=50, ftol=1e-12, verbose=False):
    if fun(a) * fun(b) > 0:
        c = None
        msg = "The function should have a sign change in the interval."
    else:
        for cont in range(niter):
            qa = fun(a)
            qb = fun(b)
            c = (a*qb - b*qa)/(qb - qa)
            qc = fun(c)
            if verbose:
                print("n: {}, c: {}".format(cont, c))
            msg = "Maximum number of iterations reached."
            if abs(qc) < ftol:
                msg = "Root found with desired accuracy."
                break
            elif qa * qc < 0:
                b = c
            elif qb * qc < 0:
                a = c
    return c, msg

def fun(x):
    return cos(x) - x

print(regula_falsi(fun, 0.5, 0.25*pi))

Julia

function regula_falsi(fun, a, b, niter=50, ftol=1e-12, verbose=false)
    if fun(a) * fun(b) > 0
        c = nothing
        msg = "The function should have a sign change in the interval."
    else
        for cont = 1:niter
            qa = fun(a)
            qb = fun(b)
            c = (a*qb - b*qa)/(qb - qa)
            qc = fun(c)
            if verbose
                println("n: $(cont), c: $(c)")
            end
            if abs(fun(c)) < ftol
                msg = "Root found with desired accuracy."
                break
            elseif qa * qc < 0
                b = c
            elseif qb * qc < 0
                a = c
            end
            msg = "Maximum number of iterations reached."
        end
    end
    return c, msg
end

function fun(x)
    return cos(x) - x
end

println(regula_falsi(fun, 0.5, 0.25*pi))

Comparison

Regarding number of lines we have: 29 in Python and 32 in Julia. The comparison in execution time is done with %timeit magic command in IPython and @benchmark in Julia.

For Python:

%timeit regula_falsi(fun, 0.5, 0.25*pi)

with result

10000 loops, best of 3: 35.1 µs per loop

For Julia:

@benchmark regula_falsi(fun, 0.5, 0.25*pi)

with result

BenchmarkTools.Trial:
  memory estimate:  48 bytes
  allocs estimate:  2
  --------------
  minimum time:     449.495 ns (0.00% GC)
  median time:      464.371 ns (0.00% GC)
  mean time:        493.785 ns (0.52% GC)
  maximum time:     9.723 μs (92.54% GC)
  --------------
  samples:          10000
  evals/sample:     198

Comments

Comments powered by Disqus