source('utils.R')

# 3.2-1:
#
root_fn = function(x){ (x-1.3333)^3 }

xs = seq( -2, +5, length.out=100 )
fs = root_fn(xs)
plot( xs, fs, type='l', xlab='x', ylab='(x-1.3333)^3' )
abline(h=0, col='black')
grid()

truth = 1.3333

mrf_res = modified_regula_falsi(1, 2, root_fn, xtol=1.e-5, ftol=1.e-16)
print(sprintf('Truth= %f; Modified Regula Falsi= %f', truth, mrf_res[1]))


# 3.2-3:
#
root_fn = function(x){ exp(x) - 1 - x - x^2/2 }

xs = seq( -0.25, +0.25, length.out=100 )
fs = root_fn(xs)
plot( xs, fs, type='l', xlab='x', ylab='(x-1.3333)^3' )
abline(h=0, col='black')
grid()


# 3.2-8:
#

# Part (a):
#
root_fn = function(x){ exp(-x) - sin(x) }
root_fn_prime = function(x){ -exp(-x) - cos(x) }

xs = seq( -1, +2, length.out=100 )
fs = root_fn(xs)
plot( xs, fs, type='l', xlab='x', ylab='exp(-X) - sin(x)')
abline(h=0, col='black')
grid()


bm_res = bisection_method(0, 1.5, root_fn, xtol=1.e-9)
mrf_res = modified_regula_falsi(0, 1.5, root_fn, xtol=1.e-9)
sec_res = secant_method(0, 1.5, root_fn, xtol=1.e-9)
new_res = newtons_method(0.5, root_fn, root_fn_prime, xtol=1.e-9)

print(sprintf('exp(-x) - sin(x)=0: %23s: x= %10f f= %10f n=%3d', 'Bisection', bm_res[1], bm_res[2], bm_res[3]))
print(sprintf('exp(-x) - sin(x)=0: %23s: x= %10f f= %10f n=%3d', 'Modified Regula Falsi', mrf_res[1], mrf_res[2], mrf_res[3]))
print(sprintf('exp(-x) - sin(x)=0: %23s: x= %10f f= %10f n=%3d', 'Secant', sec_res[1], sec_res[2], sec_res[3]))
print(sprintf('exp(-x) - sin(x)=0: %23s: x= %10f f= %10f n=%3d', 'Newton', new_res[1], new_res[2], new_res[3]))


# Part (b):
#
root_fn = function(x){ x - exp(-x^2) }
root_fn_prime = function(x){ 1 + 2*x*exp(-x^2) }

xs = seq( 0, +1.5, length.out=100 )
fs = root_fn(xs)
plot( xs, fs, type='l', xlab='x', ylab='x - exp(-x^2)')
abline(h=0, col='black')
grid()

bm_res = bisection_method(0, 1.5, root_fn, xtol=1.e-9)
mrf_res = modified_regula_falsi(0, 1.5, root_fn, xtol=1.e-9)
sec_res = secant_method(0, 1.5, root_fn, xtol=1.e-9)
new_res = newtons_method(0.5, root_fn, root_fn_prime, xtol=1.e-9)

print(sprintf('x - exp(-x^2)=0: %23s: x= %10f f= %10f n=%3d', 'Bisection', bm_res[1], bm_res[2], bm_res[3]))
print(sprintf('x - exp(-x^2)=0: %23s: x= %10f f= %10f n=%3d', 'Modified Regula Falsi', mrf_res[1], mrf_res[2], mrf_res[3]))
print(sprintf('x - exp(-x^2)=0: %23s: x= %10f f= %10f n=%3d', 'Secant', sec_res[1], sec_res[2], sec_res[3]))
print(sprintf('x - exp(-x^2)=0: %23s: x= %10f f= %10f n=%3d', 'Newton', new_res[1], new_res[2], new_res[3]))



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

xs = seq( 1, 2.0, length.out=100 )
fs = root_fn(xs)
plot( xs, fs, type='l', xlab='x', ylab='x^3 - x - 2')
abline(h=0, col='black')
grid()

bm_res = bisection_method(1.0, 1.8, root_fn, xtol=1.e-9)
mrf_res = modified_regula_falsi(1.0, 1.8, root_fn, xtol=1.e-9)
sec_res = secant_method(1.0, 1.8, root_fn, xtol=1.e-9)
new_res = newtons_method(1.5, root_fn, root_fn_prime, xtol=1.e-9)

print(sprintf('x^3 - x - 2=0: %23s: x= %10f f= %10f n=%3d', 'Bisection', bm_res[1], bm_res[2], bm_res[3]))
print(sprintf('x^3 - x - 2=0: %23s: x= %10f f= %10f n=%3d', 'Modified Regula Falsi', mrf_res[1], mrf_res[2], mrf_res[3]))
print(sprintf('x^3 - x - 2=0: %23s: x= %10f f= %10f n=%3d', 'Secant', sec_res[1], sec_res[2], sec_res[3]))
print(sprintf('x^3 - x - 2=0: %23s: x= %10f f= %10f n=%3d', 'Newton', new_res[1], new_res[2], new_res[3]))