First function is given as:
f1<-function(x){
if(x==0)
return(1)
else
return(sin(x)/x)
}Implement the trapezoidal and the Simpson’s quadrature rules to numerically compute the de nite integral above. The algorithms are presented in [3], please see Chapter 5. Hint: you can approximate the inde nite integral by considering a large interval [????a;+a] (for example a = 106). Consider equidistant nodes fxngN n=0; i.e., xn = ????a+n2a N ; n = 0; 1; : : : ;N, where N is a large integer.
TrapezoidRule = function(a, b, m, f)
{
h = (b-a)/(m-1)
Sum = 0.5 * h * (f(a)+f(b))
for (i in 1:(m-2)) {
ai = a + i*h
Sum = Sum + h*f(ai)
}
return(Sum)
}
print(TrapezoidRule(b=10^6, a=-(10^6), m=12000,f1))## [1] 3.141606SimpsonRule= function(a, b, m, f)
{
h = (b-a)/(m-1)
Sum = (1/3) * h * (f(a)+f(b))
for (i in 1:(m-2)) {
if(i%%2 != 0){
Sum = Sum + h*f(a + i*h) *(4/3)
}else{
Sum = Sum + h*f(a + i*h) *(2/3)
}
}
return(Sum)
}
print(SimpsonRule(a=-(10^6),b=10^6, m=12000,f1))## [1] 3.141626Compute the truncation error for the numerical algorithms implemented in a) for a particular a 2 R and N 2 N. That is create a function of a and N that will output IN ???? , where IN;a is the numerical approximation of the integral. Study the changes in the approximation as N and a increase as well as the di erence between the two quadrature approximations. Please write your observations.
Trapezoidal_Integration<-function(a,b,n,f)
{
table = matrix(0, nrow=n, ncol=4, dimnames=list(
c(1:n), c('m', 'Integration Result', 'Error', 'Error Ratio')))
for(i in 1:n) {
table[i,1] = ifelse(1 == i, yes=2, no=2*table[i-1,'m']-1)
table[i,2] = TrapezoidRule(a, b, table[i,'m'], f)
table[i,3] = abs(table[i,2] - pi )
table[i,4] = ifelse(1 == i, NA, table[i-1,3] /table[i,3] )
}
return(table)
}
Trapezoidal_Integration(a=-(10^6),b=10^6, n=24,f1)## m Integration Result Error Error Ratio
## 1 2 -2.099961 5.241554e+00 NA
## 2 3 999999.650006 9.999965e+05 5.241572e-06
## 3 5 500000.180666 4.999970e+05 2.000005e+00
## 4 9 249998.741156 2.499956e+05 2.000023e+00
## 5 17 125000.782187 1.249976e+05 2.000003e+00
## 6 33 62501.898891 6.249876e+04 2.000002e+00
## 7 65 31249.431008 3.124629e+04 2.000198e+00
## 8 129 15623.160000 1.562002e+04 2.000400e+00
## 9 257 7813.138512 7.809997e+03 2.000003e+00
## 10 513 3905.002313 3.901861e+03 2.001608e+00
## 11 1025 1950.932618 1.947791e+03 2.003223e+00
## 12 2049 977.035095 9.738935e+02 2.000004e+00
## 13 4097 486.947222 4.838056e+02 2.012985e+00
## 14 8193 241.903106 2.387615e+02 2.026313e+00
## 15 16385 122.522087 1.193805e+02 2.000004e+00
## 16 32769 59.690306 5.654871e+01 2.111109e+00
## 17 65537 28.274393 2.513280e+01 2.249997e+00
## 18 131073 15.707960 1.256637e+01 2.000005e+00
## 19 262145 9.424769 6.283176e+00 2.000002e+00
## 20 524289 3.141594 1.250211e-06 5.025695e+06
## 21 1048577 3.141591 1.267803e-06 9.861240e-01
## 22 2097153 3.141591 1.729308e-06 7.331270e-01
## 23 4194305 3.141591 1.837869e-06 9.409310e-01
## 24 8388609 3.141591 1.864620e-06 9.856533e-01Simpson_Integration<-function(a,b,n,f)
{
table = matrix(0, nrow=n, ncol=4, dimnames=list(
c(1:n), c('m', 'Integration Result', 'Error', 'Error Ratio')))
for(i in 1:n) {
table[i,1] = ifelse(1 == i, yes=2, no=2*table[i-1,'m']-1)
table[i,2] = SimpsonRule(a, b, table[i,'m'], f1)
table[i,3] = abs(table[i,2]- pi )
table[i,4] = ifelse(1 == i, NA, table[i-1,3] /table[i,3] )
}
return(table)
}
Simpson_Integration(a=-(10^5),b=10^5, n=20,f1)## m Integration Result Error Error Ratio
## 1 2 1.906603e-01 2.950932e+00 NA
## 2 3 1.333334e+05 1.333302e+05 2.213251e-05
## 3 5 3.333068e+04 3.332754e+04 4.000602e+00
## 4 9 1.666349e+04 1.666035e+04 2.000410e+00
## 5 17 8.334855e+03 8.331713e+03 1.999631e+00
## 6 33 4.164628e+03 4.161487e+03 2.002100e+00
## 7 65 2.085007e+03 2.081865e+03 1.998922e+00
## 8 129 1.039849e+03 1.036707e+03 2.008151e+00
## 9 257 5.225606e+02 5.194190e+02 1.995898e+00
## 10 513 2.628540e+02 2.597124e+02 1.999978e+00
## 11 1025 1.330014e+02 1.298598e+02 1.999944e+00
## 12 2049 6.387650e+01 6.073490e+01 2.138141e+00
## 13 4097 3.036803e+01 2.722644e+01 2.230733e+00
## 14 8193 1.361289e+01 1.047130e+01 2.600102e+00
## 15 16385 5.235309e+00 2.093716e+00 5.001297e+00
## 16 32769 1.046519e+00 2.095074e+00 9.993518e-01
## 17 65537 3.141820e+00 2.275482e-04 9.207164e+03
## 18 131073 3.141613e+00 2.080951e-05 1.093482e+01
## 19 262145 3.141613e+00 2.002761e-05 1.039041e+00
## 20 524289 3.141613e+00 1.998960e-05 1.001901e+00From observing the convergence characterstics, it seems like that the Simpsons method get conveged at lower values of n compared to trapezoidal
In a typical scenario we do not know the true value of the integral. Thus, to ensure the convergence of the numerical algorithms we pick a small tolerance value " and we check at every iteration k = 1; 2; : : : if the following condition holds: jIk ???? Ik????1j < “; where Ik is the value of the integral at step k. When the condition holds, the algorithm stops. Evaluate the number of steps until the algorithms from a) reach convergence for” = 10????4. What do you observe?
Trapezoidal_Integration_tol<-function(a,b,n,f,tolerance=(10^-4))
{
table = matrix(0, nrow=n, ncol=4, dimnames=list(
c(1:n), c('m', 'Integration Result', 'Error', 'Error Ratio')))
for(i in 1:n) {
table[i,1] = ifelse(1 == i, yes=2, no=2*table[i-1,'m']-1)
table[i,2] = TrapezoidRule(a, b, table[i,'m'], f)
table[i,3] = abs(table[i,2]-pi)
table[i,4] = ifelse(1 == i, NA, table[i-1,3] /table[i,3] )
if(table[i,3] < tolerance && table[i,1]>2){
print(table[i,])
print(paste("Number of Iterations=",i))
break
}
}
}#Simpson_Integration_tol<-function(a,b,n,f,tolerance=(10^-4))
# {
# table = matrix(0, nrow=n, ncol=4, dimnames=list(
# c(1:n), c('m', 'Integration Result', 'Error', 'Error Ratio')))
# for(i in 1:n) {
# table[i,1] = ifelse(1 == i, yes=2, no=2*table[i-1,'m']-1)
# table[i,2] = SimpsonRule(a, b, table[i,'m'], f1)
# table[i,3] = abs(table[i,2]- pi)
# table[i,4] = ifelse(1 == i, NA, table[i-1,3] /table[i,3] )
# if(table[i,3]<tolerance && table[i,1]>2)
# return(table[i,])
# }
# # return(tail(table,1))
#
Simpson_Integration_tol<-function(a,b,n,f,tolerance=(10^-4))
{
table = matrix(0, nrow=n, ncol=4, dimnames=list(
c(1:n), c('m', 'Integration Result', 'Error', 'Error Ratio')))
for(i in 1:n) {
table[i,1] = ifelse(1 == i, yes=2, no=2*table[i-1,'m']-1)
table[i,2] = SimpsonRule(a, b, table[i,'m'], f1)
table[i,3] = abs(table[i,2]- pi)
table[i,4] = ifelse(1 == i, NA, table[i-1,3] /table[i,3] )
if(table[i,3] < tolerance && table[i,1]>2){
print(table[i,])
print(paste("Number of Iterations=",i))
break
}
}
}Trapezoidal_Integration_tol(a=-(10^6),b=10^6, n=24,f1) ## m Integration Result Error
## 5.242890e+05 3.141594e+00 1.250211e-06
## Error Ratio
## 5.025695e+06
## [1] "Number of Iterations= 20"Simpson_Integration_tol(a=-(10^6),b=10^6, n=24,f1)## m Integration Result Error
## 1.048577e+06 3.141591e+00 2.107140e-06
## Error Ratio
## 9.939493e+05
## [1] "Number of Iterations= 21"From the above analysis we can observe that if the range of integration is large, then Trapezoidal rule is better than Simpsons, whele if it is a small range then Simpsons converge faster in minimum number of steps
f2<-function(x){
return(1+(exp(-x)*sin(8*(x^(2/3)))))
}Trapezoidal_Integration_f2<-function(a,b,n,f,tolerance=(10^-4))
{
table = matrix(0, nrow=n, ncol=4, dimnames=list(
c(1:n), c('m', 'Integration Result', 'Error', 'Error Ratio')))
for(i in 1:n) {
table[i,1] = ifelse(1 == i, yes=2, no=2*table[i-1,'m']-1)
table[i,2] = TrapezoidRule(a, b, table[i,'m'], f)
table[i,3] = abs(table[i,2]-2)
table[i,4] = ifelse(1 == i, NA, table[i-1,3] /table[i,3] )
if(table[i,3] < tolerance){
print(table[i,])
print(paste("Number of Iterations=",i))
break
}
}
}
Simpson_Integration_f2<-function(a,b,n,f,tolerance=(10^-4))
{
table = matrix(0, nrow=n, ncol=4, dimnames=list(
c(1:n), c('m', 'Integration Result', 'Error', 'Error Ratio')))
for(i in 1:n) {
table[i,1] = ifelse(1 == i, yes=2, no=2*table[i-1,'m']-1)
table[i,2] = SimpsonRule(a, b, table[i,'m'], f1)
table[i,3] = abs(table[i,2]- 2)
table[i,4] = ifelse(1 == i, NA, table[i-1,3] /table[i,3] )
if(table[i,3] < tolerance){
print(table[i,])
print(paste("Number of Iterations=",i))
break
}
}
}TrapezoidRule(a=0,b=2, m=20000,f2)## [1] 2.016279SimpsonRule(a=0,b=2, m=20000,f2)## [1] 2.016246