accounts_whole = runif(10000, 100, 1e+05)
illegal = 0
days = 0
bankillegal = function() {
while (illegal < 1e+06) {
accounts_whole = accounts_whole * (1 + 0.05/365)
accounts_trunc = floor(accounts_whole * 100)/100
illegal = (illegal * (1 + 0.05/365))
illegal = illegal + sum(accounts_whole - accounts_trunc)
days = days + 1
print(days)
}
return(days)
}
The number of days is 9630.
Given the interval [2.5, 3.5] and using the bisection method, the estimated number of interations is at least 27. It can be calculated by solving for n in the following expressions:
n = -(log(5 * 10^-9)/log(2)) - (3.5 - 2.5)
n
## [1] 26.58
Using the bisection to estimate the number of iterations:
f = function(x) {
(cos(x))^2 + sin(x) - 1
}
bisect = function(f, interval, tol, verbose = FALSE) {
a = interval[1]
b = interval[2]
fa = f(a)
fb = f(b)
j = 0 # counter for verbose printing
if (sign(fa) * sign(fb) >= 0)
stop("f(a)f(b)<0 not satisfied")
while ((b - a)/2 > tol) {
c = (a + b)/2
fc = f(c)
j = j + 1
if (verbose == TRUE) {
# show values at each step
print(j)
print(c(a, c, b, b - a))
}
if (fc == 0)
break
if (sign(fc) * sign(fa) < 0) {
b = c
fb = fc
} else {
a = c
fa = fc
}
}
return((a + b)/2)
return(j)
}
bisect(f, c(2.5, 3.5), 0.5 * 10^-8, verbose = TRUE)
## [1] 1
## [1] 2.5 3.0 3.5 1.0
## [1] 2
## [1] 3.00 3.25 3.50 0.50
## [1] 3
## [1] 3.000 3.125 3.250 0.250
## [1] 4
## [1] 3.125 3.188 3.250 0.125
## [1] 5
## [1] 3.1250 3.1562 3.1875 0.0625
## [1] 6
## [1] 3.12500 3.14062 3.15625 0.03125
## [1] 7
## [1] 3.14062 3.14844 3.15625 0.01562
## [1] 8
## [1] 3.140625 3.144531 3.148438 0.007812
## [1] 9
## [1] 3.140625 3.142578 3.144531 0.003906
## [1] 10
## [1] 3.140625 3.141602 3.142578 0.001953
## [1] 11
## [1] 3.1406250 3.1411133 3.1416016 0.0009766
## [1] 12
## [1] 3.1411133 3.1413574 3.1416016 0.0004883
## [1] 13
## [1] 3.1413574 3.1414795 3.1416016 0.0002441
## [1] 14
## [1] 3.1414795 3.1415405 3.1416016 0.0001221
## [1] 15
## [1] 3.142e+00 3.142e+00 3.142e+00 6.104e-05
## [1] 16
## [1] 3.142e+00 3.142e+00 3.142e+00 3.052e-05
## [1] 17
## [1] 3.142e+00 3.142e+00 3.142e+00 1.526e-05
## [1] 18
## [1] 3.142e+00 3.142e+00 3.142e+00 7.629e-06
## [1] 19
## [1] 3.142e+00 3.142e+00 3.142e+00 3.815e-06
## [1] 20
## [1] 3.142e+00 3.142e+00 3.142e+00 1.907e-06
## [1] 21
## [1] 3.142e+00 3.142e+00 3.142e+00 9.537e-07
## [1] 22
## [1] 3.142e+00 3.142e+00 3.142e+00 4.768e-07
## [1] 23
## [1] 3.142e+00 3.142e+00 3.142e+00 2.384e-07
## [1] 24
## [1] 3.142e+00 3.142e+00 3.142e+00 1.192e-07
## [1] 25
## [1] 3.142e+00 3.142e+00 3.142e+00 5.960e-08
## [1] 26
## [1] 3.142e+00 3.142e+00 3.142e+00 2.980e-08
## [1] 27
## [1] 3.142e+00 3.142e+00 3.142e+00 1.490e-08
## [1] 3.142
The general function of the Newton method for this function is as follow: \[ x_{i+1}=x_{i}-f(x_{i})/f'{x_{i+1}} \]
\[ x_{i+1}=x_{i}-[(cos(x_{i}))^2+sin(x_{i})-1]/[-2cos(x_{i})sin(x_{i})+cos(x_{i})] \]
The convergence rate to \( \pi \) with a starting guess of 3.1 is quadratic.
x_0 = 3.1
vec <- rep(0, 50)
NMf = function(x) {
return(x - ((cos(x))^2 + sin(x) - 1)/(-2 * cos(x) * sin(x) + cos(x)))
}
for (i in 1:50) {
x_0 = NMf(x_0)
vec[i] = x_0
}
abserror = abs(vec - pi)
errorval <- rep(0, 49)
errorval2 <- rep(0, 49)
for (i in 2:length(abserror)) {
errorval[i] = abserror[i]/abserror[i - 1]
errorval2[i] = abserror[i]/(abserror[i - 1])^2
}
plot(errorval[2:50])
plot(errorval2[2:50])
Repeating d, but with starting guess of 1.5. Yes, the answer is different from d. The convergence rate to \( \pi/2 \) is linear, as \( e_{i+1}/e_{i} \) converges to a value between 0 and 1. For d, the convergence rate to \( \pi \) is quadratic, as \( e_{i+1}/e_{i}^2 \) is positive.
x_0 = 1.5
vec <- rep(0, 50)
NMf = function(x) {
return(x - ((cos(x))^2 + sin(x) - 1)/(-2 * cos(x) * sin(x) + cos(x)))
}
for (i in 1:50) {
x_0 = NMf(x_0)
vec[i] = x_0
}
abserror = abs(vec - pi/2)
errorval <- rep(0, 49)
errorval2 <- rep(0, 49)
for (i in 2:length(abserror)) {
errorval[i] = abserror[i]/abserror[i - 1]
errorval2[i] = abserror[i]/(abserror[i - 1])^2
}
plot(errorval[2:50], type = "l")
plot(errorval2[2:50], type = "l")
f = function(x) {
det(rbind(c(1, 2, 3, x), c(4, 5, x, 6), c(7, x, 8, 9), c(x, 10, 11, 12))) -
1000
}
options(digits = 6)
bisect(f, c(7, 11), 0.5 * 10^-6, verbose = FALSE)
## [1] 9.7083
The value of x is 9.7083.
library(Matrix)
## Loading required package: lattice
f = function(x) {
H = Hilbert(5)
H[1, 1] = x
return(max(eigen(H)$values) - pi)
}
options(digits = 6)
bisect(f, c(0, 10), 0.5 * 10^-6, verbose = FALSE)
## [1] 2.94801