An example
fy = function(x)
{
5*x^6 - 36*x^5 + 165/2*x^4 - 60*x^3 + 36
}
fdy = function(x)
{
30*x^5 - 180*x^4 + 330*x^3 - 180*x^2
}
fd2y = function(x)
{
150*x^4 - 720*x^3 + 990*x^2 - 360*x
}
plot(fy, 0, 3.3, ylim=c(-50,50))
x = seq(0, 3.3, 0.1)
y = fdy(x) / 2
lines(x, y, lty=2)
abline(h=0)
next.x = function(x)
{
x - fdy(x) / fd2y(x)
}
max.iter =15
tol = 0.001
xi = matrix(nrow=max.iter, ncol=5)
xi[1,2] = 0.7
for(i in 1:length(xi[,1])) {
xi[i,1] = i
xi[i,3] = fdy(xi[i,2])
xi[i,4] = fd2y(xi[i,2])
xi[i,5] = next.x(xi[i,2])
if(abs(xi[i,2] - xi[i,5]) < tol) break
else xi[i+1,2] = xi[i,5]
}
xi
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0.7000 -1.319e+01 22.16 1.2952
## [2,] 2 1.2952 1.785e+01 52.24 0.9535
## [3,] 3 0.9535 -2.717e+00 56.64 1.0015
## [4,] 4 1.0015 8.750e-02 60.09 1.0000
## [5,] 5 1.0000 6.306e-05 60.00 1.0000
## [6,] NA NA NA NA NA
## [7,] NA NA NA NA NA
## [8,] NA NA NA NA NA
## [9,] NA NA NA NA NA
## [10,] NA NA NA NA NA
## [11,] NA NA NA NA NA
## [12,] NA NA NA NA NA
## [13,] NA NA NA NA NA
## [14,] NA NA NA NA NA
## [15,] NA NA NA NA NA
tol = 0.000000001
xi = matrix(nrow=max.iter, ncol=4)
xi[1,1] = 1.6
for(i in 1:length(xi[,1])) {
xi[i,2] = fdy(xi[i,1])
xi[i,3] = fd2y(xi[i,1])
xi[i,4] = next.x(xi[i,1])
if(abs(xi[i,1] - xi[i,4]) < tol) break
else xi[i+1, 1] = xi[i,4]
}
xi
## [,1] [,2] [,3] [,4]
## [1,] 1.600 2.580e+01 -7.68 4.960
## [2,] 4.960 1.696e+04 25498.56 4.295
## [3,] 4.295 5.420e+03 10714.92 3.789
## [4,] 3.789 1.697e+03 4601.51 3.421
## [5,] 3.421 5.075e+02 2070.60 3.175
## [6,] 3.175 1.357e+02 1036.75 3.045
## [7,] 3.045 2.645e+01 649.51 3.004
## [8,] 3.004 2.078e+00 548.97 3.000
## [9,] 3.000 1.685e-02 540.07 3.000
## [10,] 3.000 1.139e-06 540.00 3.000
## [11,] 3.000 -9.095e-13 540.00 3.000
## [12,] NA NA NA NA
## [13,] NA NA NA NA
## [14,] NA NA NA NA
## [15,] NA NA NA NA