# 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]
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!
# 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]
$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} $$# 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]