math4610
Tasksheet 5
Author: David Merkley
Language: Python
Task 1
Routine Name: NewtonRoot
Description/Purpose: Approximates a root with Newton’s method. I added this to my python package.
Input: A function f(x), the derivative of f(x), f’(x), starting point x0, a tolerance, and a max iteration.
Output: Approximation of the root
Implementation/Code:
from math import exp
def f(x):
return x * exp(3 * x**2) - 7 * x
def derivf(x):
return ((6 * x**2) + 1) * exp(3 * x**2) - 7
def newtonMethod(f, derivf, x0, tol=1e-10, maxiter=1e5):
fx = f(x0)
error = 10 * tol
counter = 0
while counter in range(int(maxiter)) and error > tol:
x = x0 - fx / derivf(x0)
error = abs(x - x0)
fx = f(x)
if error < tol:
return x
x0 = x
print("No root was found")
print(newtonMethod(f, derivf, 1))
Output:
0.8053798584219568
Task 2
Routine Name: SecantRoot
Description/Purpose: Approximates a root with the Secant method. I added this to my python package.
Input: A function f(x), starting point x0, initial point x1a tolerance, and a max iteration.
Output: Approximation of the root
Implementation/Code:
from math import exp
def f(x):
return x * exp(3 * x**2) - 7 * x
def derivf(x):
return ((6 * x**2) + 1) * exp(3 * x**2) - 7
def secantMethod(f, x0, x1, tol=1e-10, maxiter=1e10):
xk1 = x1
xk0 = x0
f1 = f(xk1)
f0 = f(xk0)
error = 10 * tol
counter = 0
while counter in range(int(maxiter)) and error > tol:
error = abs(xk1 - xk0)
xk2 = xk1 - f1 / ((f1 - f0) / (xk1 - xk0))
xk0 = xk1
xk1 = xk2
f0 = f1
f1 = f(xk1)
if error < tol:
return xk1
print("No root was found")
print(secantMethod(f, 1, 2))
Output:
0.8053798584219568
Task 3
Routine Name: NewtonConvergence
Description/Purpose: Verifies quadratic convergence using Newton’s method.
Input: A function f(x), the derivative of f(x), f’(x), starting point x0, a tolerance, and a max iteration.
Output: A log-log plot
Implementation/Code:
import math
from math import exp
from matplotlib import pyplot as plt
def f(x):
return x * exp(3 * x**2) - 7 * x
def derivf(x):
return ((6 * x**2) + 1) * exp(3 * x**2) - 7
def logplot(logerr):
x = []
y = []
for i in range(1, len(logerr)):
y.append(logerr[i])
x.append(logerr[i - 1])
plt.title("Log-Log Plot of Newton's Convergence")
plt.xlabel("Log $e_k$")
plt.ylabel("log $e_{k - 1}$")
plt.plot(x, y)
plt.grid()
plt.show()
def newtonConvergence(f, derivf, x0, tol=1e-10, maxiter=1e10):
x = [x0]
fx0 = f(x[0])
counter = 1
while counter in range(1, int(maxiter + 1)):
error = abs(fx0)
tempx = (x[counter - 1] - fx0 / derivf(x[counter - 1]))
fx0 = f(tempx)
if error < tol:
logerr = [math.log(abs(xcounter - tempx)) for xcounter in x]
logplot(logerr)
return
x.append(tempx)
counter = counter + 1
newtonConvergence(f, derivf, 1)
Output:
Task 4
Routine Name: SecantConvergence
Description/Purpose: Verifies quadratic convergence using the Secant method.
Input: A function f(x), the derivative of f(x), f’(x), starting point x0, a tolerance, and a max iteration.
Output: A log-log plot
Implementation/Code:
import math
from math import exp
from matplotlib import pyplot as plt
def f(x):
return x * exp(3 * x**2) - 7 * x
def logplot(logerr):
x = []
y = []
for i in range(1, len(logerr)):
y.append(logerr[i])
x.append(logerr[i - 1])
plt.title("Log-Log Plot of Secant's Convergence")
plt.xlabel("Log $e_k$")
plt.ylabel("log $e_{k - 1}$")
plt.plot(x, y)
plt.grid()
plt.show()
def secantConvergence(f, x0, x1, tol=1e-10, maxiter=1e10):
x = [x0, x1]
fx0 = f(x[0])
fx1 = f(x[1])
counter = 2
while counter in range(2, int(maxiter + 2)):
# error = abs(fx0)
tempx = x[counter - 1] - fx0 / ((fx0 - fx1) / (x[counter - 1] - x[counter - 2]))
fx1 = fx0
fx0 = f(tempx)
if abs(fx0) < tol:
logerr = [math.log(abs(xcounter - tempx)) for xcounter in x]
logplot(logerr)
return
x.append(tempx)
counter = counter + 1
secantConvergence(f, 1, 2)
Output:
Task 5
Routine Name: HybridMethod
Description/Purpose: Searches for roots by combining bisection method with the Newton’s method.
Input: A function f(x), the derivative of f(x), f’(x), starting point x0, a tolerance, and a max iteration.
Output: A log-log plot
Implementation/Code:
import math
from math import exp
from matplotlib import pyplot as plt
def f(x):
return x * exp(3 * x**2) - 7 * x
def derivf(x):
return ((6 * x**2) + 1) * exp(3 * x**2) - 7
def newton(f, fx, derivf, x):
return x - fx / derivf(x)
def bisection(f, fa, fb, a, b):
# Different signs for f(a) and f(b)
if fa > 0 and fb < 0 or fa < 0 and fb > 0:
counter = 0
while counter in range(5):
c = (a + b) / 2
fc = f(c)
if (fa * fc < 0):
fb = fc
b = c
else:
fa = fc
a = c
counter = counter + 1
c = (a + b) / 2
return a, b, c
else:
print("Values were not opposite signs")
def hybridMethod(f, derivf, a, b, tol=1e-10, maxiter=1e10):
c = (a + b) / 2
fa = f(a)
fb = f(b)
fc = f(c)
for counter in range(int(maxiter)):
c = newton(f, fc, derivf, c)
if c < a or c > b:
a, b, c = bisection(f, fa, fb, a, b)
fc = f(c)
else:
fc = f(c)
if fa * fc < 0:
fb = fc
b = c
else:
fa = fc
a = c
error = b - a
if error < tol:
return c
print("No root found")
print(hybridMethod(f, derivf, -5, 10))
Output:
0.8053798584219568
Task 6
The Bisection method is a root finding method that is very slow. It’s simple and usually used to obtain a rough approximation of a root. Newton’s method has some problems, because you need to calcuate the derivative of the function before even starting to use the method. Though this method is extremely fast. The secant method has a low number of steps, but is slower than Newton’s method. Also, the secant method does not require that the root remain bracketed.