Integration by Parts

When find the solution through the Simpson rule for "I" we must have to choose the step size neither small nor large to get the accurate answer. It means we have to limit the number of sub segments for the higher accuracy. If we select large number of sub segments it should affect to increase the round off error, and if we choose small number of sub segments we can’t get exact solution.

There I used multiple segments Simpson 1/3 rule and Simpson 3/8 rule to solve above equation. Because Simpson 1/3 only works for even number of segments and Simpson 3/8 works for multiplication of 3 segments.

Multiple segments Simpson 1/3 rule

Multiple segments Simpson 3/8 rule

According to my program user can select the number of segments (n) as they want. When we apply a = 0, b = 1, assume n = 50 segments, h = (1-0 )/50 we can do the calculation in below progress. Here is the result from Simpson 1/3 rule.

 

 When above constant values apply to Simpson 3/8 rule we can get the solution which is very closer to Simpson 1/3 rule. When compared the solutions of the recursion method is little bit of far away from the Simpson result. Here is the recursion scenario. There n should be 10.

Using recursive formula evaluates the solution through the FORTRAN program it does not give the answer which comes from Simpson rules.

We can calculate absolute relative true error as following.

Implementation in FORTRAN

Non Recursion Solution

1.)Do the compilation as below.

    gfortran -ffree-form simpsonint_portion.f95
    ./a.out

2.)User must have to input the lower limit and upper limit of the integrations.

    E.g.: Lower limit : 0
   
          Upper limit : 1

3.)Select the choice Simpson1/3(enter-1) or Simpson3/8(enter-2)

    E.g.: Enter the choice : 1

4.)Enter the number of sub intervals.

    E.g.: Enter the number of sub intervals : 50

5.)Now you can see the outputs on a Terminal.

    E.g.: Answer of I(n) from the Simpson1/3 rule is :    0.08388

simpsonint_portion.f95 Code

program simpsonint_portion

Integer :: interval, i, choice
real :: a, b, h, rslt = 0.0
real :: evensum = 0.0, oddsum = 0.0, sum = 0.0
!Take the inputs from the user to identify the integration range
print *, ''
print *, "Enter the lower limit and upper limit of the integrations....................."
write(*,2,advance="no") "Lower limit : "
read *, a
write(*,2,advance="no") "Upper limit : "
read *, b

print *,"------------------------Choose an option---------------------------------"
print *, '1. Simpson 1/3 rule'
print *, '2. Simpson 3/8 rule'
print *,"-------------------------------------------------------------------------"

print *, ''

!Take the choice from the user
write (*,2,advance="no") "Enter the choice : "
2 format('',A20)
read *,choice

print *, ''

!Take the number of sub segments from the user
write (*,3,advance="no") "Enter the number of sub intervals : "
3 format('',A37)
read *, interval
!Calculation of the step size(h)
h = (b-a)/interval

if(choice == 1 .and. mod(interval,2) == 0)then
    rslt = simp13(a,b)    !Call to the simpson 1/3 rd rule
    print *,"-------------------------------------------------------------------------"
    write(*,10)rslt
    10 format(1x,'Answer of I(n) from the Simpson1/3 rule is : ',F10.5)
    print *,"-------------------------------------------------------------------------"
else if(choice == 2 .and. mod(interval,3) == 0)then
    rslt = simp38(a,b)    !Call to the simpson 3/8 rule
    print *,"-------------------------------------------------------------------------"
    write(*,20)rslt
    20 format(1x,'Answer of I(n) from the simpson3/8 rule is : ',F10.5)
    print *,"-------------------------------------------------------------------------"
else
    print *, "Invalid input.................."
    stop
end if

contains

!Simpson 1/3 rd function
function simp13(a,b)
 
    real :: simp13, a, b
    do i=1, (interval-1)
    if(mod(i,2) == 0)then
      evensum = evensum + f((a+(i*h)))
    else
      oddsum = oddsum + f((a+(i*h)))
    end if
    end do
    !Simpson 1/3 rd formula
    simp13 = (h/3) * (f(a) + 2*evensum + 4*oddsum + f(b))

end function simp13

!Simpson 3/8 function 
function simp38(a,b)
 
    real :: simp38, a, b
    do i=1, (interval-1)
    if(mod(i,3) == 0)then
      sum = sum + 2*f((a+(i*h)))
    else
      sum = sum + 3*f((a+(i*h)))
    end if
    end do
    !Simpson 3/8 formula
    simp38 = (3*h/8) * (f(a) + sum + f(b))
 
end function simp38

end program simpsonint_portion

!Function of f(x)
function f(x)

    implicit none
    real :: f, x, y=1.0
    f = x**10 * exp(x-y)
 
end function f

Recursion Solution

1.)Do the compilation as below.

    gfortran -ffree-form recursion.f95
    ./a.out

2.)User must have to input the power of a function.

    E.g.: Enter the power(n) of the function : 10

3.)Now you can see the outputs on a Terminal.

    E.g.: Answer of recursion of I(n) is :    0.05067

recursion.f95 Code

program recursion

    real :: rslt, x=1.0
    Integer :: n
    !Take an input from the user
    print *,'.........................................................'
    write (*,2,advance="no") "Enter the power(n) of the function : "
    2 format('',A38)
    read *, n
    !Call to the function solve
    rslt = solve(n)
    !print *, rslt
    write(*,20)rslt
    20 format(1x,'Answer of recursion of I(n) is : ',F10.5)
    print *,'.........................................................'

contains

!Recursive function to calculate I(n)
recursive function solve(n) result(a)
 
    Integer :: n
    real :: a
    if(n == 1)then
      a = (1.0/exp(x))
      return
    else
      a = (1 - real(n)*solve(n-1)) !Recursion portion
      return
    end if

end function solve
 
end program recursion