First Steps in Julia: Iterations, Newton's method, Eigensolver and Shank's transform¶

In [1]:
# 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]

roots and eigenvalues:¶

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 [2]:
# 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])

[1.0000000000000029, 1.9999999999999971, 2.9999999999999885, 4.000000000000014]
[-1.865174681370263e-14, -7.105427357601002e-15, 0.0, 1.1368683772161603e-13]

Shanks transformation¶

$S_N$ denotes a series summed up to N terms, then

$$ S = \frac{S_{N+1}S_{N-1}-S_N^2}{S_{N+1}+S_{N-1}−2 S_N} $$
In [3]:
# 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]
In [ ]: