Computational Methods for Numerical Analysis with R
Numerical Analysis Computation with R is an overview of numerical analysis topics using R. This shows how common topics in numerical analysis such as interpolation, numerical integration, roots of non linear equation (using bisection method, newton raphson method), finite difference, newton forward and backward difference, and the likes can be implemented in R. Every algorithm described is given with complete function implementation in R, along with examples to demostrate the function and its use. This project is working with the assumption that we all have basic knowledge of topics in numerical analysis, and presume we are not new to R. So, let's start with solution to non-linear equation using bisection method.
1.0 Bisection Method
The bisection method is used to find the root of a countinuous function \(f(x)\) on an interval \([a,b]\). The principle behind this method is the intermidiate theorem for continuous functions. It works by narrowing the gap between positive and negative intervals until it closes to the correct answer. This method narrows the gap by taking the average of the positive and negative intervals.
If \(a\) and \(b\) are such that \(f(a)\) and \(f(b)\) have opposite signs, there is at least one root of the equation \(f(x) = 0\) in the interval between \(a\) and \(b\) (provided \(f(x)\) is a continuous function). In the method of bisection, the midpoint of the interval \(c = \frac{a + b}{2}\) is taken, and from the sign of \(f(c)\) it can be deduced whether a root lies in the half interval to the left or right of \(c\). Whichever half interval is indicated, its midpoint is then taken and the procedure repeated. This method often requires many iterations, therefore it can be rather slow but never fails to produce the root. The procedure stops when two values of x are equal to the required degree of accuracy.
1.1 The Bisection Method in R
In R, there are two ways of finding the roots of a countinuous function using Bisection Method. One of which is a function in NLRoot package in R and the other is by writing the algorithm for the bisection method.
To use the BFfzero function in R, we have to install the package NLRoot first, then read the library as follows.
library(NLRoot)Example 1
Suppose we are to find the root(s) of the function \(x^3-10x+4 = 0\). It is advisable to first plot the function to know where the root(s) possibly lie(s).
func = function(x){
x^3 - 10*x + 4
}
curve(func, xlim=c(-4,4), col = 'red', lwd = 2, lty = 2, xlab = "x", ylab = "f(x)")
abline(h=0)
abline(v=0)
It is obvious from the curve that there are three roots for the function. To find these roots, we use BFfzero function below in R.
BFfzero(func, -4, -2)## [1] 1
## [1] -3.345963
## [1] 1.819209e-05
## [1] "finding root is successful"BFfzero(func, 0, 2)## [1] 1
## [1] 0.4067291
## [1] -6.637231e-06
## [1] "finding root is successful"BFfzero(func, 2, 3)## [1] 1
## [1] 2.939236
## [1] 2.518246e-05
## [1] "finding root is successful"These roots are -3.345963, 0.4067291 and 2.939236
Example 2
We wish to find the root of the function \(sinx+x-1 = 0\) in this example, We plot the function to find an interval containing the root(s), as follows:
func2 = function(x){
sin(x) + x - 1
}
curve(func2, xlim=c(-4,4), col = 'red', lwd = 2, lty=2, xlab = "x", ylab = "f(x)")
abline(h=0)
abline(v=0)
Then, we find the root by using the interval \([0,2]\)
library(NLRoot)
BFfzero(func2, 0, 2)## [1] 1
## [1] 0.5109711
## [1] -4.418654e-06
## [1] "finding root is successful"Now, we have used the R inbuilt function for the bisection method in finding the root of non-linear equation. Though, there are inbuilt functions in R for other numerical computation such as Newton-Raphson Method, Difference Equation, Simpson rule for integration and the likes, our interest is to write a code to define a function for these numerical computations. Below is the one for bisection method.
Bisection <- function(f, a, b, n = 1000, tol = 1e-7) {
if (!(f(a) < 0) && (f(b) > 0)) {
stop('The root does not exist within this interval')
} else if (!(f(a) > 0) && (f(b) < 0)) {
stop('The root does not exist within this interval')
}
for (i in 1:n) {
c <- (a + b) / 2
if ((f(c) == 0) || ((b - a) / 2) < tol) {
return(c)
}
ifelse(sign(f(c)) == sign(f(a)),
a <- c,
b <- c)
}
print('Too many iterations')
}Now, the function Bisection() is ready to be used, but before that, there is need to define the function to find its root as follows.
f = function(x){
x^3 - 10*x + 4
}After defining the function as above,we use the function Bisection(f, a, b) to find its roots. Where a and b are the points on the curve and f is the defined function \(f(x)\).
Bisection(f, -4, -2)## [1] -3.345963Bisection(f, 0, 2)## [1] 0.4067284Bisection(f, 2, 3)## [1] 2.939235Using this function Bisection(), these roots obtained for \(f(x)\) are the same as that of the inbuilt function BFfzero.
With this, we can use Bisection(f, a, b) to find the root of a defined function in R. Where a and b are the points on the curve, where there is a change in sign and f is the defined function \(f(x)\).
2.0 Newton-Raphson Iterative Method
Newton-Raphson iterative method is another way of approximating the root of non-linear continuous functions. It requires the derivative of the function in question and the initial value \(x_0\).
Consider the graph \(y = f(x)\), the \(x\) - value at point A when the graph crosses the \(x\) - axis is given as the solution of the equation \(f(x) = 0\). If P is the point near to the curve, then \(x = x_0\) is an approximated value.
If the intial value (\(x_0\)) of \(x\) is known, the value of \(x_1\) can be determined, and this process can be repeated till the root converges at a point. The iterative procedure is defined as \[x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\]
Suppose we are to find the root of the equation \(x^2 - 5x - 40 = 0\), when \(x_0 = 2\), the code and the output are given below.
nrm = function(f, df, x0){
n = 0
while(n < 100){
n = n+1
y0 = f(x0)
y1 = df(x0)
if(y1 < 0.0001){
break(0)
}
x1 = x0 - (y0/y1)
if(abs((x0-x1)/x1)<0.0001){
print("The point converges here")
break(0)
}
x0 = round(x1,4)
print(paste("Iteration" ,n, "is",x0))
}
}The code above defined the function nrm(f, df, x), where f, df and x are \(f(x)\), \(f'(x)\) and initial value \(x_0\) respectively. Now, let's define the function \(f(x)\) and its derivative \(f'(x)\) as follows.
f = function(x){
z = x^3 - 5*x - 40
return(z)
}
df = function(x){
z = 3*x^2 - 5
return(z)
}nrm means Newton-Raphson Method. To find the root of a defined function using Newton-Ramphson Method in R, use nrm(f, df, x)
nrm(f, df, 2)## [1] "Iteration 1 is 8"
## [1] "Iteration 2 is 5.6898"
## [1] "Iteration 3 is 4.4333"
## [1] "Iteration 4 is 3.9706"
## [1] "Iteration 5 is 3.9057"
## [1] "Iteration 6 is 3.9044"
## [1] "The point converges here"The point of convergence is our root, hence the root of the equation at \(x_0 = 4\) is 3.9044 to 4 d.p. To find other roots of the equation, we use another intial condition \(x_0\), based on the point where there is a change in sign.
3.0 Lagrange Interpolation
Interpolation may be described as the art of reading between the lines in a table. We are concerned with obtaining values of a function at some points when the values of the function is known at some other points.
The lagrange interpolating formular is given as:
\[P_n(x) = \sum\limits_{k=0}^n \frac{L_k(x)}{L_k(x_k)}.f(x_k)\] Where the lagrange basis functions is given as:
\[\frac{L_k(x)}{L_k(x_k)} = \frac{(x - x_0)(x - x_1)...(x-x_j)(x-x_{k-1})}{(x_j - x_0)(x_j - x_1)...(x_j-x_{j-1})(x_j-x_{j-1})}\]
Here is how we defined a function for Lagrange Interpolation using the interpolating formular above, 'the commented code' here gives a prompt that allow users to input the values to be computed. To make this happen, uncomment the code and source the code, it will prompt the user to input their own values.
#n = as.integer(readline(prompt = "Enter n: "))
#for(k in 1:n)
#{
# x = as.double(readline(prompt = "Enter x: "))
# y = as.double(readline(prompt = "Enter y: "))
# }
#a = as.integer(readline(prompt = "Enter the value of a: "))
f = 0.0
Lagrange = function(x, y, n, a){
for (i in 1:n){
LX = 1.0
LXk = 1.0
for (j in 1:n){
if(i != j){
LX = LX*(a - x[j])
LXk = LXk*(x[i] - x[j])
}
j = j+1
}
f = f + (LX/LXk)*y[i]
i = i+1
}
s = round(f, 4)
print(paste("f(",a,") is ", s))
}Suppose \(f(x) = \ln(x)\), and we wish to estimate the value of \(\ln(0.6)\) from the table below.
| x | 0.4000 | 0.5000 | 0.7000 | 0.8000 |
| y | -0.9165 | -0.6931 | -0.3567 | -0.2231 |
You assign the variables as follows;
x = c(0.4, 0.5, 0.7, 0.8)
y = c(-0.9165, -0.6931, -0.3567, -0.2231)But, if the commented codes are 'uncommented', there is no need to assign variables, all you need is a command prompt that allows you to input your variables one after the other.
Now, after assigning the variables, use the function Lagrange(x, y, n, a)
Lagrange(x, y, 4, 0.66)## [1] "f( 0.66 ) is -0.4151"The output here gives the estimated value of \(\ln(0.66)\) to be -0.4151 using the given table for lagrange interpolation. Comparing the result with the actual value of \(\ln(0.66)\) which is -0.4155, gives a relative error of 0.096.
The interpolating values can be plotted using:
plot(x,y, col = "green", lwd = 1, pch = 19, xlab = "x", ylab = "f(x)")
yb = approx(x,y, xout = 0.66)
points(0.66, yb$y, col = "red", lwd = 1, pch = 19)
4.0 Newton Forward Difference
Newton forward difference formular is used to interpolate the values of \(y\) that is closer to the beginning of a set of tabular values. \(y_0\) may be taken as any point of the table, but the formular contains only those values of \(y\) which comes after the value chosen as \(y_0\) The Newton forward difference formular is given as: \[f(x) = y_0 + r\Delta y_0 + \frac{r(r-1)}{2!}\Delta^2y_0 + \frac{r(r-1)(r-2)}{3!}\Delta^3y_0 + \frac{r(r-1)(r-2)(r-3)}{4!}\Delta^4y_0 + ...\] Where \[x = x_0 + rh\] and \[r = \frac{x - x_0}{h}\]
#n = as.integer(readline(prompt = "Enter n: "))
##x = c(rep(0,n))
#y = c(rep(0,n))
#for(i in 1:n)
#{
# x[i] = as.double(readline(prompt = "Enter x: "))
# y[i] = as.double(readline(prompt = "Enter y: "))
#}
nfd = function(x,y,a){
n = length(x)
A = matrix(c(rep(0,(n-1)^2)), nrow = n-1, ncol = n-1, byrow = TRUE)
for (j in 1:(n-1)){
for (i in 1:(n-j))
{
if(j==1)
A[i,j]=(y[i+1]-y[i])
else
A[i,j]=(A[(i+1),(j-1)]-A[i,(j-1)])
}
}
#no = as.integer(readline(prompt = "Enter number to find: "))
u = (a - x[1])/(x[2]-x[1])
print(A)
sum = y[1]
for (i in 1:(n-1))
{
prod = 1
for(j in 0:(i-1))
prod = prod*(u-j)
prod = prod/factorial(i)
sum = sum + prod*A[1,i]
}
print(sum)
}If the population of a town was given in the table below and we are to estimate the population for year 1895 using Newton forward difference formular,
| Year | 1891 | 1901 | 1911 | 1921 | 1931 |
| Population | 46 | 66 | 81 | 93 | 101 |
then we define the table by assigning Year to x and Population to y as a vector, as shown below.
x = c(1891,1901,1911,1921,1931)
y = c(46,66,81,93,101)
a = 1895The function nfd(x, y, a) is then used to estimate the population, using the three parameters, where a is 1895.
nfd(x,y,a)## [,1] [,2] [,3] [,4]
## [1,] 20 -5 2 -3
## [2,] 15 -3 -1 0
## [3,] 12 -4 0 0
## [4,] 8 0 0 0
## [1] 54.85285.0 Newton Backward Difference
Newton's forward difference formular cannot be used for interpolating a value of y near the end of table of values. For this purpose, we use another formular known as Newton backward difference formular. and it can be written as: \[f(x) = y_n + r\nabla y_n + \frac{r(r+1)}{2!}\nabla^2y_n + \frac{r(r+1)(r+2)}{3!}\nabla^3y_n + \frac{r(r+1)(r+2)(r+3)}{4!}\nabla^4y_n + ...\] Where
\[x = x_0 + rh\]
and \[r = \frac{x - x_0}{h}\]
#n = as.integer(readline(prompt = "Enter n: "))
#x = c(rep(0,n))
#y = c(rep(0,n))
#for(ii in 1:n)
#{
# x[ii] = as.double(readline(prompt = "Enter x: "))
# y[ii] = as.double(readline(prompt = "Enter y: "))
#}
nbd = function(x,y,a){
n = length(x)
A = matrix(c(rep(0,(n)^2)), nrow = n, ncol = n, byrow = TRUE)
for (j in 1:(n-1)){
for (i in (j+1):n)
{
if(j==1)
A[i,j]=(y[i]-y[i-1])
else
A[i,j]=(A[(i),(j-1)]-A[(i-1),(j-1)])
}
}
#no = as.integer(readline(prompt = "Enter no to find: "))
u = (a - x[n])/(x[2]-x[1])
print(A)
sum = y[n]
for (i in 1:(n-1))
{
prod = 1
for(j in 0:(i-1))
prod = prod*(u+j)
prod = prod/factorial(i)
sum = sum + prod*A[n,i]
}
print(sum)
}Using same tabular values used for the Newton forward difference, estimate for the year 1925 using Newton backward difference.
| x | 1891 | 1901 | 1911 | 1921 | 1931 |
| y | 46 | 66 | 81 | 93 | 101 |
Following the same procedure as the Newton forward difference, we have:
x = c(1891,1901,1911,1921,1931)
y = c(46,66,81,93,101)
a = 1925nbd(x,y,a)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0 0 0
## [2,] 20 0 0 0 0
## [3,] 15 -5 0 0 0
## [4,] 12 -3 2 0 0
## [5,] 8 -4 -1 -3 0
## [1] 96.83686.0 Trapezoidal Rule for Integration
Trapezoidal rule is an integration rule that approximates the area under the curve by dividing the total area into smaller trapezoids, then find the area of the small trapezoids in a definite interval. This rule is used for approximating the definite integrals where linear approximations of the function are used. There are other methods of approximating definite integrals, such as Simpson Rule, Rectanglar rule and the likes, which are going to be discussed later. To calculate the area under the curve using the general trapezoidal rule, we use: \[\int_a^bf(x)dx = \frac{h}{2}[(y_0 + y_n) + 2(y_1 + y_2 + y_3 + ... + y_{n-1})]\] This trapezoidal general formular helps in defining a function for trapezoidal rule approximation for integration.
#a = as.integer(readline(prompt = "Enter first value: "))
#b = as.integer(readline(prompt = "Enter last value: "))
#n = as.integer(readline(prompt = "Enter the interval n: "))
trap = function(f,a,b,n){
h = (b-a)/n
x = seq(a, b, by = h)
f = sapply(x, f)
t = h*(f[1]/2 + sum(f[2:n]) + f[n+1]/2)
A = signif(t,4)
print(paste("The integral of f(x) on interval [",a ,"," ,b, "] when n =" ,n ,"is" ,A))
}Suppose we are to approximate the integral \(\int_1^7 \frac{x^3}{1+x^4}dx\) using n = 30, the function \(f(x) = \frac{x^3}{1+x^4}dx\) has to be defined first. To define the function, we use:
f = function(x){
x^3*(1+x^4)^-1
}The function trap(f,a,b,n) can then be used to approximate the integral of the defined function as follows.
trap(f,1,7,30)## [1] "The integral of f(x) on interval [ 1 , 7 ] when n = 30 is 1.771"The bigger the n, the greater the degree of accuracy to the actual integral.
7.0 Simpson's Rule for Integration
Simpson's Rule is another rule for approximating definite integrals, it converges faster compared to the Trapezoid rule.
Simpson's rule is defined as:
\[\int_a^bf(x)dx = \frac{h}{3}[(y_0 + y_n) + 2(y_2 + y_4 + y_6 + ... + y_{n-2}) + 4(y_1 + y_3 + y_5 + ... + y_{n-1})]\] Where \[h = \frac{b-a}{n} \] We can write the function to approximate a definite integral using Simpson's rule as follows:
#a = as.integer(readline(prompt = "Input first value: "))
#b = as.integer(readline(prompt = "Input last value: "))
#n = as.integer(readline(prompt = "Input the interval n: "))
simp = function(f,a,b,n){
x[1]=f(a)
x[n+1]=f(b)
h = (b-a)/n
s = f(a) + f(b)
for(i in 1:(n-1))
{
x[i-1]=f(a+i*h)
if(i%%2!=0){
s = s + 4*f(a+i*h)
}
else{
s = s + 2*f(a+i*h)
}
}
I = signif(s*(h/3),4)
print(paste("The integral of f(x) on interval [",a ,"," ,b, "] when n =" ,n ,"is" ,I))
}Using the same function we integrated using trapezoid rule, we have:
f = function(x){
x^3*(1+x^4)^-1
}simp(f,1,7,30)## [1] "The integral of f(x) on interval [ 1 , 7 ] when n = 30 is 1.773"Summary
R as an integrated suite that facilitates data manipulation, calculation, mathematical computation and graphical display. All these have been used as a tool to compute some numerical analysis problems. Algorithms were defined to write functions to compute each numerical method discussed early with examples to implement the functions in R.
References
G.Shanker Rao (2003). Numerical Analysis (3rd ed.). New Age International Publishers.
Simpson's rule for integration (2017). In RPubs. Retrieved from https://rpubs.com/aaronsc32/simpsons-rule-approximate-integrals
Bisection method (2016). In Wikipedia. Retrieved from https://en.wikipedia.org/wiki/Bisection_Method