source('utils.R')

# Exercise 3.5-3:
#

# Part (a):
#
root_fn = function(x){ exp(-x) - x }
root_fn_p = function(x){ -exp(-x) - 1 }

# Plot the function (to see approximately where the root is):
#
xs = seq(0, 1, length.out=100)
ys = root_fn(xs)
plot(xs, ys, type='l')
grid()

#
# Decide on an interval surrounding the root and check that the conditions of Theorem 3.2 are satisfied
#
a = 0
b = 1
print( root_fn(a) * root_fn(b) ) # negative
xs = seq(a, b, length.out=100)
yps = root_fn_p(xs)
plot(xs, yps, type='l') # f'(x) is never zero and f"(x) > 0 on this interval
print( c( abs(root_fn(a)/root_fn_p(a)), abs(root_fn(b)/root_fn_p(b)), b-a ) ) # both less than b-a = 1
n_res = newtons_method(0.5*(a+b), root_fn, root_fn_p)
s_res = secant_method(a, b, root_fn)
print(sprintf('Part (a): Newton x_0= %.3f; x= %f f= %f n= %d', 0.5*(a+b), n_res[1], n_res[2], n_res[3]))
print(sprintf('Part (a): Secant [a, b]= [%.1f, %.1f]; x= %f f= %f n= %d', a, b, s_res[1], s_res[2], s_res[3]))


# Part (b):
#
root_fn = function(x){ x^3 - x - 1 }
root_fn_p = function(x){ 3*x^2 - 1 }

# Plot the function (to see approximately where the root is):
#
xs = seq(1, 2, length.out=100)
ys = root_fn(xs)
plot(xs, ys, type='l')
grid()

# Decide on an interval surrounding the root and check that the conditions of Theorem 3.2 are satisfied
#
a = 1
b = 2
print( root_fn(a) * root_fn(b) ) # negative
xs = seq(a, b, length.out=100)
yps = root_fn_p(xs)
plot(xs, yps, type='l') # f'(x) is never zero and f"(x) > 0 on this interval
print( c( abs(root_fn(a)/root_fn_p(a)), abs(root_fn(b)/root_fn_p(b)), b-a ) ) # both less than b-a = 1
n_res = newtons_method(0.5*(a+b), root_fn, root_fn_p)
s_res = secant_method(a, b, root_fn)
print(sprintf('Part (b): Newton x_0= %.3f; x= %f f= %f n= %d', 0.5*(a+b), n_res[1], n_res[2], n_res[3]))
print(sprintf('Part (b): Secant [a, b]= [%.1f, %.1f]; x= %f f= %f n= %d', a, b, s_res[1], s_res[2], s_res[3]))


# Part (c):
#
root_fn = function(x){ exp(-x^2) - cos(x) }
root_fn_p = function(x){ - 2 * x * exp(-x^2) + sin(x) }
root_fn_pp = function(x){ - 2 * exp(-x^2) + 4 * x^2 * exp(-x^2) + cos(x) }

# Plot the function (to see approximately where the root is):
#
xs = seq(1, 2, length.out=100)
ys = root_fn(xs)
plot(xs, ys, type='l')
grid()

# Decide on an interval surrounding the root and check that the conditions of Theorem 3.2 are satisfied
#
a = 1.2
b = 1.6
xs = seq(a, b, length.out=100)
print( root_fn(a) * root_fn(b) ) # negative
yps = root_fn_p(xs)
plot(xs, yps, type='l') # f'(x) is never zero
ypps = root_fn_pp(xs)
plot(xs, ypps, type='l') # f"(x) > 0 on this interval
print( c( abs(root_fn(a)/root_fn_p(a)), abs(root_fn(b)/root_fn_p(b)), b-a ) ) # both less than b-a
n_res = newtons_method(0.5*(a+b), root_fn, root_fn_p)
s_res = secant_method(a, b, root_fn)
print(sprintf('Part (c): Newton x_0= %.3f; x= %f f= %f n= %d', 0.5*(a+b), n_res[1], n_res[2], n_res[3]))
print(sprintf('Part (c): Secant [a, b]= [%.1f, %.1f]; x= %f f= %f n= %d', a, b, s_res[1], s_res[2], s_res[3]))


# Exercise 3.5-6:
#
root_fn = function(x){ x - tan(x) }
root_fn_p = function(x){ 1 - 1/cos(x)^2 }
root_fn_pp = function(x){ 2*sin(x)/cos(x)^3 }

# Plot the function (to see approximately where the root is):
#
a = 102.09
b = 102.0999
xs = seq(a, b, length.out=100)
ys = root_fn(xs)
plot(xs, ys, type='l', xlab='x', ylab='y', main='y = x - tan(x)')
grid ()


# Decide on an interval surrounding the root and check that the conditions of Theorem 3.2 are satisfied
#
print( root_fn(a) * root_fn(b) ) # negative
yps = root_fn_p(xs)
plot(xs, yps, type='l', xlab='x', ylab='y prime', main='first derivative') # f'(x) is never zero
ypps = root_fn_pp(xs)
plot(xs, yps, type='l', xlab='x', ylab='y prime prime', main='second derivative') # and f"(x) < 0 on this interval
print( c( abs(root_fn(a)/root_fn_p(a)), abs(root_fn(b)/root_fn_p(b)), b-a ) ) # both less than b-a
n_res = newtons_method(0.5*(a+b), root_fn, root_fn_p)
s_res = secant_method(a, b, root_fn)
print(sprintf('Ex 3.5-6: Newton x_0= %.3f; x= %f f= %f n= %d', 0.5*(a+b), n_res[1], n_res[2], n_res[3]))
print(sprintf('Ex 3.5-6: Secant [a, b]= [%.1f, %.1f]; x= %f f= %f n= %d', a, b, s_res[1], s_res[2], s_res[3]))