This notebook contains an excerpt from the Python Programming And Numerical Methods: A Guide For Engineers And Scientists; the content is available on GitHub.

The text is released under the CC-BY-NC-ND license, and code is released under the MIT license. If you find this content useful, please consider supporting the work by buying the book!

< 19.4 Newton-Raphson Method | Contents | CHAPTER 20. Numerical Differentiation >

Summary¶

Roots are an important property of functions.
The bisection method is a way of finding roots based on divide and conquer. Although stable, it might converge slowly compared to the Newton-Raphson method.
The Newton-Raphson method is a different way of finding roots based on approximation of the function. The Newton-Raphson method converges quickly close to the actual root, but can have unstable behavior.

Problems¶

Write a function \(my\_nth\_root(x, n, tol)\), where \(x\) and \(tol\) are strictly positive scalars, and \(n\) is an integer strictly greater than 1. The output argument, \(r\), should be an approximation \(r = \sqrt[N]{x}\), the \(N\)-th root of \(x\). This approximation should be computed by using the Newton Raphson method to find the root of the function \(f(y) = y^N - x\). The error metric should be \(|f(y)|\).
Write a function \(my\_fixed\_point(f, g, tol, max_iter)\), where \(f\) and \(g\) are function objects and \(tol\) and \(max\_iter\) are strictly positive scalars. The input argument, \(max\_iter\), is also an integer. The output argument, \(X\), should be a scalar satisfying \(|f(X) - g(X)| < tol\); that is, \(X\) is a point that (almost) satisfies \(f(X) = g(X)\). To find \(X\), you should use the Bisection method with error metric, \(|F(m)| < tol\). The function \(my\_fixed\_point\) should “give up” after \(max\_iter\) number of iterations and return \(X = []\) if this occurs.
Why does the bisection method fail for \(f(x) = 1/x\) with error given by \(|b-a|\)? Hint: How does \(f(x)\) violate the intermediate value theorem?
Write a function \(my\_bisection(f, a, b, tol)\), that returns \([R, E]\), where \(f\) is a function object, \(a\) and \(b\) are scalars such that \(a < b\), and \(tol\) is a strictly positive scalar value. The function should return an array, \(R\), where \(R[i]\) is the estimation of the root of \(f\) defined by \((a + b)/2\) for the \(i\)-th iteration of the bisection method. Remember to include the initial estimate. The function should also return an array, \(E\), where \(E[i]\) is the value of \(| f(R[i]) |\) for the \(i\)-th iteration of the bisection method. The function should terminate when \(E(i) < tol\). You may assume that \({\text{sign}}(f(a)) \neq {\text{sign}}(f(b))\).

Clarification: The input \(a\) and \(b\) constitute the first iteration of bisection, and therefore \(R\) and \(E\) should never be empty.

Test Cases:

In: f = lambda x: x**2 - 2
    [R, E] = my_bisection(f, 0, 2, 1e-1)
Out: R = [1, 1.5, 1.25, 1.375, 1.4375]
     E = [1, 0.25, 0.4375, 0.109375, 0.06640625]
       
In: f = lambda x: np.sin(x) - np.cos(x)
    [R, E] = my_bisection(f, 0, 2, 1e-2)
Out: R = [1, 0.5, 0.75, 0.875, 0.8125, 0.78125]
     E = [0.30116867893975674, 0.39815702328616975, 0.05005010885048666, 0.12654664407270177, 0.038323093040207645, 0.005866372111545948]

Write a function \(my\_newton(f, df, x0, tol)\), that returns \([R, E]\), where \(f\) is a function object, \(df\) is a function object to the derivative of \(f\), \(x0\) is an initial estimation of the root, and \(tol\) is a strictly positive scalar. The function should return an array, \(R\), where \(R[i)]\) is the Newton-Raphson estimation of the root of \(f\) for the \(i\)-th iteration. Remember to include the initial estimate. The function should also return an array, \(E\), where \(E[i]\) is the value of \(| f(R[i]) |\) for the \(i\)-th iteration of the Newton-Raphson method. The function should terminate when \(E(i) < tol\). You may assume that the derivative of \(f\) will not hit 0 during any iteration for any of the test cases given.

Test Cases:

In: f = lambda x: x**2 - 2
   df = lambda x: 2*x
    [R, E] = my_newton(f, df, 1, 1e-5)
Out: R = [1, 1.5, 1.4166666666666667, 1.4142156862745099]
     E = [1, 0.25, 0.006944444444444642, 6.007304882871267e-06]
       
In: f = lambda x: np.sin(x) - np.cos(x)
   df = lambda x: np.cos(x) + np.sin(x)
    [R, E] = my_newton(f, df, 1, 1e-5)
Out: R = [1, 0.782041901539138, 0.7853981759997019]
     E = [0.30116867893975674, 0.004746462127804163, 1.7822277875723103e-08]

Consider the problem of building a pipeline from an offshore oil platform, a distance \(H\) miles from the shoreline, to an oil refinery station on land, a distance \(L\) miles along the shore. The cost of building the pipe is \(C_{\text{ocean}/\text{mile}}\) while the pipe is under the ocean and \(C_{\text{land}/\text{mile}}\) while the pipe is on land. The pipe will be built in a straight line toward the shore where it will make contact at some point, \(x\), between 0 and \(L\). It will continue along the shore on land until it reaches the oil refinery. See the figure for clarification.

Write a function \(my\_pipe\_builder(C\_ocean, C\_land, L, H)\), where the input arguments are as described earlier, and \(x\) is the x-value that minimizes the total cost of the pipeline. You should use the bisection method to determine this value to within a tolerance of \(1\times10^{-6}\) starting at an initial bound of \(a=0\) and \(b=L\).

Test Cases:

In: my_pipe_builder(20, 10, 100, 50)
Out: 28.867512941360474
  
In: my_pipe_builder(30, 10, 100, 50)
Out: 17.677670717239380
  
In: my_pipe_builder(30, 10, 100, 20)
Out: 7.071067392826080

Find a function \(f(x)\) and guess for the root of \(f, x_0\), such that the Newton-Raphson method would oscillate between \(x_0\) and \(-x_0\) indefinitely.

< 19.4 Newton-Raphson Method | Contents | CHAPTER 20. Numerical Differentiation >

Python Numerical Methods

Summary¶

Problems¶