Numerical Analysis

bisection.jl

function plot_bisection_method(f, a, b;
                               n_max_iter=20, n_samples=100, interval=1000)
    xrange, yrange, xs, ys = init_xs_ys(f, a, b, n_samples)
    fig, ax = init_figure_plot()

    function draw(i)
        cla()

        ax[:set_xlabel]("x")
        ax[:set_ylabel]("f(x)")
        ax[:set_xlim](xrange)
        ax[:plot](xs, zeros(length(xs)), color="black", alpha=0.6)  # x-axis
        ax[:plot](xs, ys, color="blue")

        c = (a + b) / 2

        ax[:vlines]([a], min(0, f(a)), max(0, f(a)), color="black", alpha=0.6)
        ax[:vlines]([b], min(0, f(b)), max(0, f(b)), color="black", alpha=0.6)
        ax[:scatter]([c], [0], color="red")

        if f(a)*f(c) < 0
            b = c
        else
            a = c
        end
    end

    @pyimport matplotlib.animation as animation
    animation.FuncAnimation(fig, draw,
                            frames=1:n_max_iter,
                            interval=interval)
end


function bisection_method(f, a, b; n_max_iter = 200, threshold = 0.01)
    """
    Find ξ such that `f(ξ) = 0` in the interval `[a, b]` using Bisection method
    """
    if f(a)*f(b) >= 0
        throw("Invalid initial values (a = $a b = $b)")
    end

    xs = Float64[]
    for i in 1:n_max_iter
        c = (a + b) / 2

        push!(xs, c)
        if abs(f(c)) < threshold
            return xs
        end

        if f(a)*f(c) < 0
            b = c
        else
            a = c
        end
    end
    return xs
end

comparison.jl

ENV["PYTHON"] = rstrip(readstring(`which python3`))

using Formatting
using ForwardDiff
using PyCall
using PyPlot


include("plotlib.jl")
include("bisection.jl")
include("newton.jl")
include("secant.jl")
include("false_position.jl")

f(x) = x.^2 - 2

a, b = 0, 5

n_max_iter = 10
x = b  # inital value for Newton method
writer = "imagemagick"

anim = plot_bisection_method(f, a, b; n_max_iter = n_max_iter)
anim[:save]("bisection.gif", writer = writer)

anim = plot_newton_method(f, x, a, b; n_max_iter = n_max_iter)
anim[:save]("newton.gif", writer = writer)

anim = plot_secant_method(f, a, b; n_max_iter = n_max_iter)
anim[:save]("secant.gif", writer = writer)

anim = plot_false_position_method(f, a, b; n_max_iter = n_max_iter)
anim[:save]("false_position.gif", writer = writer)


xs_bisection = bisection_method(f, a, b)
xs_newton = newton_method(f, x)
xs_secant = secant_method(f, a, b)
xs_false_position = false_position_method(f, a, b)

target = sqrt(2)
xs_bisection = abs.(xs_bisection - target)
xs_newton = abs.(xs_newton - target)
xs_secant = abs.(xs_secant - target)
xs_false_position = abs.(xs_false_position - target)

N = max(
    length(xs_bisection),
    length(xs_newton),
    length(xs_secant),
    length(xs_false_position)
)
ts = 1:N


fig, ax = init_figure_plot()
ax[:plot](1:length(xs_bisection), xs_bisection,
          label = "Bisection", linestyle = "solid")
ax[:plot](1:length(xs_newton), xs_newton,
          label = "Newton", linestyle="dashed")
ax[:plot](1:length(xs_secant), xs_secant,
          label = "Secant", linestyle="dotted")
ax[:plot](1:length(xs_false_position), xs_false_position,
          label = "False-Position", linestyle="dashdot")
ax[:set_ylabel]("Absolute difference of \$x\$ and the solution")
ax[:set_xlabel]("Time step")
ax[:legend]()
fig[:savefig]("convergence.png")

false_position.jl

function plot_false_position_method(f, a, b;
                                    n_max_iter=20, n_samples=100, interval=1000)
    xrange, yrange, xs, ys = init_xs_ys(f, a, b, n_samples)
    fig, ax = init_figure_plot()

    function draw(i)
        cla()

        ax[:set_xlabel]("x")
        ax[:set_ylabel]("f(x)")
        ax[:set_xlim](xrange)
        ax[:plot](xs, zeros(length(xs)), color="black", alpha=0.6)  # x-axis
        ax[:plot](xs, ys, color="blue")
        ax[:plot]([a, b], [f(a), f(b)], color="red")
        ax[:vlines]([a], min(0, f(a)), max(0, f(a)), color="black", alpha=0.6)
        ax[:vlines]([b], min(0, f(b)), max(0, f(b)), color="black", alpha=0.6)

        c = b - f(b) * (b-a) / (f(b)-f(a))

        if f(a) * f(c) > 0
            a = c
        else
            b = c
        end

        ax[:scatter]([c], [0], color="red")
    end

    @pyimport matplotlib.animation as animation
    animation.FuncAnimation(fig, draw,
                            frames=1:n_max_iter,
                            interval=interval)
end


function false_position_method(f, a, b; n_max_iter = 200, threshold = 0.01)
    if f(a)*f(b) >= 0
        throw("Invalid initial values (a = $a b = $b)")
    end

    xs = Float64[]
    for i in 1:n_max_iter
        c = b - f(b) * (b-a) / (f(b)-f(a))

        push!(xs, c)
        if abs(f(c)) < threshold
            return xs
        end

        if f(a) * f(c) > 0
            a = c
        else
            b = c
        end
    end
    return xs
end

jacobi.jl

function jacobi(A, b, n_iter=100)
    G = A .* (ones(size(A)...) - eye(size(A)...))
    D = A .* eye(size(A)...)

    x = rand(size(b))
    for i in 1:n_iter
        x = inv(D) * (b - G*x)
        println("x = $x")
    end
    return x
end


A = [
    1 1;
    2 5;
]
b = [3, 9]
x = jacobi(A, b)
println("A = $A")
println("b = $b")
println("x = $x")
println("A*x = $(A*x)")

newton.jl

function plot_newton_method(f, x, a, b;
                            n_max_iter=20, n_samples=100, interval=1000)
    xrange, yrange, xs, ys = init_xs_ys(f, a, b, n_samples)
    fig, ax = init_figure_plot()

    function draw(i)
        cla()

        ax[:set_xlabel]("x")
        ax[:set_ylabel]("f(x)")
        ax[:set_xlim](xrange)
        ax[:plot](xs, zeros(length(xs)), color="black", alpha=0.6)  # x-axis
        ax[:plot](xs, ys, color="blue")

        d = ForwardDiff.derivative(f, x)  # g = df / dx evaluated at x
        xₙ = x - f(x) / d

        ax[:vlines]([x], min(0, f(x)), max(0, f(x)), color="black", alpha=0.6)
        ax[:vlines]([xₙ], min(0, f(xₙ)), max(0, f(xₙ)), color="black", alpha=0.6)
        ax[:plot]([x, xₙ], [f(x), 0], color="red", linestyle="--", alpha=0.6)
        ax[:scatter]([xₙ], [0], color="black", alpha=0.6)

        x = xₙ
    end

    @pyimport matplotlib.animation as animation
    animation.FuncAnimation(fig, draw,
                            frames=1:n_max_iter,
                            interval=interval)
end


function newton_method(f, x; n_max_iter = 200, threshold = 0.01)
    """
    Find ξ such that `f(ξ) = 0` with initial value `x` using Newton method
    """
    xs = Float64[]
    for i in 1:n_max_iter
        d = ForwardDiff.derivative(f, x)  # g = df / dx evaluated at x
        x = x - f(x) / d

        push!(xs, x)
        if abs(f(x)) < threshold
            return xs
        end
    end
    return xs
end

plotlib.jl

function init_xs_ys(f, a, b, N; margin = 0.1)
    xrange = [a - margin, b + margin]
    xs = linspace(xrange..., N)
    ys = f(xs)
    yrange = [minimum(ys), maximum(ys)]
    return xrange, yrange, xs, ys
end


function init_figure_plot(figsize = (10, 10))
    fig = figure(figsize = figsize)
    ax = fig[:add_subplot](1, 1, 1)
    fig, ax
end

secant.jl

function plot_secant_method(f, a, b;
                            n_max_iter=20, n_samples=100, interval=1000)
    xrange, yrange, xs, ys = init_xs_ys(f, a, b, n_samples)
    fig, ax = init_figure_plot()

    x, xₚ = a, b
    function draw(i)
        cla()

        ax[:set_xlabel]("x")
        ax[:set_ylabel]("f(x)")
        ax[:set_xlim](xrange)
        ax[:plot](xs, zeros(length(xs)), color="black", alpha=0.6)  # x-axis
        ax[:plot](xs, ys, color="blue")

        xₙ = x - f(x) * (x-xₚ) / (f(x)-f(xₚ))

        L, H = min(xₚ, x, xₙ), max(xₚ, x, xₙ)
        ax[:vlines]([xₚ], min(0, f(xₚ)), max(0, f(xₚ)), color="black", alpha=0.6)
        ax[:vlines]([xₙ], min(0, f(xₙ)), max(0, f(xₙ)), color="black", alpha=0.6)
        ax[:vlines]([x], min(0, f(x)), max(0, f(x)), color="red")
        ax[:plot]([L, H], [f(L), f(H)], color="red", linestyle="--", alpha=0.6)

        xₚ = x
        x = xₙ
    end

    @pyimport matplotlib.animation as animation
    animation.FuncAnimation(fig, draw,
                            frames=1:n_max_iter,
                            interval=interval)
end


function secant_method(f, a, b; n_max_iter = 200, threshold = 0.01)
    """
    Find ξ such that `f(ξ) = 0` in the interval `[a, b]` using Secant method
    """
    xs = Float64[]
    x, xₚ = a, b

    for i in 1:n_max_iter
        x, xₚ = x - f(x) * (x-xₚ) / (f(x)-f(xₚ)), x
        push!(xs, x)
        if abs(f(x)) < threshold
            return xs
        end
    end
    return xs
end