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 numerical analysis topics 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 guild assumed we all have a basic knowledge of topics in numerical analysis, and presumed we are not new to R. So, let 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 interval 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 to the left or right of \(c\). Whichever half interval is indicated its midpoint is then taken and procedure repeated. The method often require many iterations, and is therefore slow but never fails to eventually 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 the 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"The roots here 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 thus.
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've 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 too, but our interest is to write a code to define a function for these numerical computations. The one below is 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, one have to define the function to find it's root as follows.
f = function(x){
x^3 - 10*x + 4
}After defining the function as above, then use the function Bisection(f, a, b) to find it's roots. Where a and b are the point 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(), the roots obtained for \(f(x)\) is 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 point 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 we are to find it's root 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 determine, 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 intial value \(x_0\) respectively. Now, Let define the function \(f(x)\) and it 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 the other roots of the equation, we use another intial condition \(x_0\), base on the point where there's a change in sign.
3.0 Lagrange Interpolation
Interpolation may be describe as the art of reading between the lines in a table. We are concern with obtaining values of a function at some points when the values of that 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 happens, 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 uncomment, there is no need to assign variables, all you need is a command prompt that allow 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\) near 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 come 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 we 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 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 approximate 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 it uses the linear approximations of the function. There are other method of approximating definte integral, such as Simpson Rule, Rectanglar rule and the likes, which are going to be discussed later here. 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 define 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 use 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 more the 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 compare 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 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 facilities data manipulation, calculation, mathematical computation and graphical display has been used as a tool to compute some numerical analysis problems. Algorithm was defined to write a function to compute each numerical methods 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