1. What is Newton-Method?
Since early education in Mathematics, we have been exposed to technique of finding roots. For instance when we need to find which value of \(y\) that can fulfills the following equation
$$ y^2 = 25,$$
we can use root finding method to solve it. Analytically, we can then proceed to solve the above equation as follows
$$ y^2 - 25 = 0$$
$$ y^2 - 5^2 = 0$$
$$ (y-5)(y+5) = 0.$$
Thus, the \(y\)'s that can become the part of our solution are 5 and -5.
This simple procedure is a hard problem when the degree of polynomial increases or involving complex functions such as trigonometric function. To make it worse, the difficulties magnified when we want to set/automate the root finding procedure into the computer. Thus, the Newton-Method which is among the earliest algorithm suggested by Joseph Raphson to approximate the solution of the roots can be used.
2. How to perform Newton Method?
We will first look at how this algorithm is derived. Recalled that in calculus, the derivative of a function \(f(x)\) is defined as,
$$ f'(x) = \lim_{h \to 0} \frac{f(x+h)-f(x)}{(x+h)-x}.$$
Our goal is to develop an algorithm that will proposed a new \(x\) in each iteration and the value of \(f(x) \rightarrow 0 \) (a root) as the iteration continues. But how do we generate the new \(x\)?
This is obvious if we manipulate the definition above as follows.
Let \(x_{n+1} = x+h\) be the new value that we want to propose and \(x_n = x\) be the former value that we already know its \(f(x_n)\) value. Then it is noted that,
$$ f'(x_n) \approx \frac{f(x_{n+1}) - f(x_n)}{x_{n+1} - x_n}.$$
$$ x_{n+1} - x_n \approx \frac{f(x_{n+1}) - f(x_n)}{f'(x_n)}$$
$$ x_{n+1} \approx x_n + \frac{f(x_{n+1}) - f(x_n)}{f'(x_n)}.$$
Since we are expecting the new \(x_{n+1}\) to cause \(f(x_{n+1})\) approaching 0, thus the last equation is updated into,
$$ x_{n+1} \approx x_n + \frac{0 - f(x_n)}{f'(x_n)} = x_n - \frac{f(x_n)}{f'(x_n)}.$$
Thus completing our algorithm.
EG:
a) Find a positive \(y\) that satisfy \(y^2 = 25\).
Solution: Denote that \(y^2 - 25 = 0\). We can use the root finding algorithm discussed above. Information that we need are as follows,
$$ f(y) = y^2 - 25$$
$$ f'(y) = 2y.$$
R script to run Newton-Raphson method:
f <- function(y) y^2 - 25
f.prime <- function(y) 2*y
y <- 3 # set initial y (close to answer is better)
tol <- 100 # initialize threshold to exit loop when f(y) is 0
n.iter <- 0 # store iteration
while(abs(tol) > 1e-5 && n.iter < 500){
tol <- f(y)
y <- y - f(y)/f.prime(y)
n.iter = n.iter + 1
}
print(paste0("With ", n.iter, " iterations, y = ",y))
## [1] "With 5 iterations, y = 5"
f <- function(x) 2*x^3 - 4*x - 1
rootSolve::uniroot.all(f, c(-5,5))
## [1] -1.2669017 -0.2586498 1.5256917
No comments:
Post a Comment