Basic Numerical Methods¶

Newton, Halley and Steffensen¶

In [169]:
# finding square root

a = pi
f = x -> x^2 - a

ret1 = 10
ret2 = 10
ret3 = 10

for i1 in 1:20
    
    ret1 = 0.5 * (ret1 + a/ret1)  # Newton's
    ret2 = (ret2^3 + 3*ret2*a)/(3*ret2^2+a)  # Halley's
    ret3 = ret3 - f(ret3)^2/(f(ret3+f(ret3)) - f(ret3)) # Steffensen
    
    println([i1, ret1, ret2, ret3])
    
end

[1.0, 5.15707963267949, 3.6096919925406867, 9.171147292301468]
[2.0, 2.883130095853957, 1.9193000519367729, 8.355843451113838]
[3.0, 1.9863883041004537, 1.772676962522595, 7.556246904888214]
[4.0, 1.7839742445187232, 1.7724538509063996, 6.775053652960133]
[5.0, 1.7724910486031316, 1.772453850905516, 6.01568834295974]
[6.0, 1.7724538512958334, 1.7724538509055159, 5.282587921744467]
[7.0, 1.7724538509055159, 1.772453850905516, 4.581636313835007]
[8.0, 1.7724538509055159, 1.7724538509055159, 3.920852559575508]
[9.0, 1.7724538509055159, 1.772453850905516, 3.311508069895036]
[10.0, 1.7724538509055159, 1.7724538509055159, 2.76992730713633]
[11.0, 1.7724538509055159, 1.772453850905516, 2.320020351575515]
[12.0, 1.7724538509055159, 1.7724538509055159, 1.9943524965935278]
[13.0, 1.7724538509055159, 1.772453850905516, 1.8211034897592875]
[14.0, 1.7724538509055159, 1.7724538509055159, 1.7753021300796064]
[15.0, 1.7724538509055159, 1.772453850905516, 1.7724642124798888]
[16.0, 1.7724538509055159, 1.7724538509055159, 1.7724538510431627]
[17.0, 1.7724538509055159, 1.772453850905516, 1.7724538509055159]
[18.0, 1.7724538509055159, 1.7724538509055159, 1.772453850905516]
[19.0, 1.7724538509055159, 1.772453850905516, 1.7724538509055159]
[20.0, 1.7724538509055159, 1.7724538509055159, 1.772453850905516]

root finding and eigenvalue solver¶

Finding eigenvalues can be done by solving for roots of characteristic polynomial. Can one go about finding all roots by using an eigenvalue solver? YES!

In [170]:
# finding all roots of a polynomial
using LinearAlgebra
using Plots

# x^n + c(n-1) x^(n-1) + ... + c(0)
# collect cvec = [c(0), ..., c(n-1)]

# roots: [1, 2, 3, 4]
cvec = [24, -50, 35, -10]

cvec = rand(3)

ndim = size(cvec)[1]
f = x -> sum([cvec[i1]*x^(i1-1) for i1 in 1:ndim]) + x^ndim

# construct the  Frobenius companion
cmat = ones(ndim, ndim)

for i1 in 1:ndim, i2 in 1:ndim
    cmat[i1, i2] = i1-i2==1 ? 1 : 0
end
cmat[:, ndim] = -cvec[:]

rvec = eigen(cmat).values

println(rvec)
println([f(rvec[i1]) for i1 in 1:ndim])

ComplexF64[-1.0295105551416812 + 0.0im, 0.1328706860625224 - 0.9564876555460818im, 0.1328706860625224 + 0.9564876555460818im]
ComplexF64[2.220446049250313e-16 + 0.0im, -7.216449660063518e-16 - 4.440892098500626e-16im, -7.216449660063518e-16 + 4.440892098500626e-16im]

Shanks transformation¶

$S_N$ means summing a series up to N terms, then

$ S( = \frac{S_{N+1}S_{N-1}-S_N^2}{S_{N+1}+S_{N-1}-2S_{N}} $

In [171]:
# Shanks

f = nn -> 1/nn^2
s1 = N -> sum([f(nn) for nn in 1:N])

function shanks(_s1, N)
    tmp1, tmp2, tmp3 = _s1(N+1), _s1(N), _s1(N-1)
    ret = (tmp1*tmp3-tmp2^2) / (tmp1+tmp3-2*tmp2)
    return ret
end

ref = pi^2/6

for Nmax in 20:120
    
    mod(Nmax,10)==0 ? 
    println([Nmax, 
            s1(Nmax), 
            shanks(s1, Nmax), 
            shanks(nn1->shanks(s1, nn1), Nmax),
            shanks(nn2->shanks(nn1->shanks(s1, nn1), nn2), Nmax),
            ref
            ]) : nothing
    
end

[20.0, 1.5961632439130233, 1.6205534878167467, 1.632748617859053, 1.6388464364800461, 1.6449340668482264]
[30.0, 1.6121501176015975, 1.6285435602232496, 1.6367402818181183, 1.6408395574635124, 1.6449340668482264]
[40.0, 1.6202439630069354, 1.632589642006862, 1.63876249317556, 1.6418327862643305, 1.6449340668482264]
[50.0, 1.6251327336215293, 1.6350337237471557, 1.639984198365324, 1.6424868790199454, 1.6449340668482264]
[60.0, 1.6284055175083685, 1.636669980252294, 1.6408023451637441, 1.6425515607902512, 1.6449340668482264]
[70.0, 1.6307499074900194, 1.6378421060799178, 1.6413882004325655, 1.642992576546274, 1.6449340668482264]
[80.0, 1.6325118663375642, 1.6387230464418254, 1.6418286360058412, 1.6434802600065364, 1.6449340668482264]
[90.0, 1.6338844555141359, 1.6394093174488094, 1.6421721233576085, 1.6425890172823152, 1.6449340668482264]
[100.0, 1.6349839001848927, 1.6399590248063813, 1.6424448636945808, 1.642338600536995, 1.6449340668482264]
[110.0, 1.635884354854303, 1.6404092416968246, 1.6426717267563733, 1.6417601209676917, 1.6449340668482264]
[120.0, 1.6366353592878375, 1.640784736756058, 1.6428617967700516, 1.6428893662404283, 1.6449340668482264]
In [172]:
# Shanks

f = nn -> 4*(-1)^nn/(2*nn+1)
s1 = N -> sum([f(nn) for nn in 0:N])


function shanks(_s1, N)
    tmp1, tmp2, tmp3 = _s1(N-1), _s1(N), _s1(N+1)
    ret = (tmp1*tmp3-tmp2^2) / (tmp1+tmp3-2*tmp2)
    return ret
end

ref = pi

for Nmax in 20:200
    
    mod(Nmax,20)==0 ? 
    println([Nmax, 
            s1(Nmax), 
            shanks(s1, Nmax), 
            shanks(nn1->shanks(s1, nn1), Nmax),
            shanks(nn2->shanks(nn1->shanks(s1, nn1), nn2), Nmax),
            ref
            ]) : nothing
    
    
end

[20.0, 3.189184782277595, 3.141565734658557, 3.1415926990259124, 3.141592661442751, 3.141592653589793]
[40.0, 3.1659792728432166, 3.1415890289407717, 3.141592655226769, 3.1415923087266373, 3.141592653589793]
[60.0, 3.157984995168666, 3.141591552545737, 3.141592653993129, 3.141592479790599, 3.141592653589793]
[80.0, 3.1539378622726155, 3.141592183260289, 3.141592652605248, 3.141592934495151, 3.141592653589793]
[100.0, 3.151493401070991, 3.1415924109720126, 3.1415926548692172, 3.1415925739966877, 3.141592653589793]
[120.0, 3.149856975293274, 3.141592512483335, 3.1415926567679504, 3.141592621049056, 3.141592653589793]
[140.0, 3.148684762993838, 3.141592564412247, 3.141592659547041, 3.1415925652326995, 3.141592653589793]
[160.0, 3.147803773812001, 3.1415925936877973, 3.1415926593666517, 3.141592652368605, 3.141592653589793]
[180.0, 3.147117473309497, 3.1415926114310815, 3.1415926625804063, 3.141592641147515, 3.141592653589793]
[200.0, 3.1465677471829565, 3.1415926228048625, 3.1415926412470925, 3.141592641483803, 3.141592653589793]