More Program Examples

Numerical Computation

Numerical Computation

•     first application of computers

•     evaluate definite integrals

•     solve differential equations

•     sum of series

•     find digits of constants such as pi

•     useful in many engineering and scientific applications involving complex mathematical  models

Numerical Integration

•     evaluate complex definite integrals

•     find area under a curve

•     curve specified by a complicated function f

Trapezoidal Method

•     split the interval of integration into many small intervals

•     in each small interval [p,q] approximate the curve by a straight line segment joining (p,f(p)) and (q,f(q))

•     add up all the areas of the resulting trapezoids ( taking into account the sign)

•     gives an approximation to area of curve

•     lengths of intervals need not be same

•     errors in approximation

•     large value of 2nd derivative causes errors

•     use simplest form

–   read number of intervals (n)

–   divide interval of integration [a,b] into intervals of length (b-a)/real(n) = h

–   endpoints of ith interval are a+(i-1)h and a+ih

•     integral = h(f(a)/2.0 + f(a+h) + f(a+2h) + …

                  + f(a+(n-1)h) +f(b)/2.0)      

Integration Program

program integrate

! integrate sqrt(1-x*x) from –1.0 to +1.0

implicit none

real :: h, integral

integer :: n, i

print *, “type number of intervals”

read *, n

h = 2.0/real(n)

integral = 0.0       ! f(-1.0) = f(1.0) = 0.0

do i = 1,n-1

integral = integral +                            &

                sqrt(1.0-(real(i)*h-1.0)**2)

end do

integral = integral * h

print *, “integral is”, integral

end program integrate

 

Output of Integration Program

•      1.5707963268.. is the correct value ( pi/2 )

type number of intervals

10

integral is   1.5185245

type number of intervals

100

integral is   1.5691345

type number of intervals

1000

integral is   1.5707436

type number of intervals

10000

integral is   1.5707954

type number of intervals

100000

integral is   1.5708289

type number of intervals

1000000

integral is   1.5707991

Results

•     two types of errors

–   truncation error due to trapezoid approximation

–   round-off error due to finite precision

•     truncation error is large if number of intervals is too small

•     round-off error is large if it is too large

•     n should be chosen carefully

Series Summation

•     many complicated functions can be written as a convergent Taylor series

•     functions can be approximately evaluated by taking finite number of terms from the series

•     the Taylor series can also be integrated or differentiated under certain conditions

•     series must be convergent for the range of values considered

Example of Summing Series

•     consider the function  sqrt(1-x2), -1 < x < 1

•     can be expanded using binomial theorem

•     f(x) = 1 – a0 – a1 – a2 – a3 - ……

–   a0 = x2/2.0

–   an = an-1* x2 *(2n-1)/(2n+2), n >= 1

•     express an in terms of an-1 rather than an explicit expression

•     series converges for |x| < 1

Program to Sum Series

program sum_series

! evaluates sqrt(1-x*x) using binomial theorem

implicit none

real :: sum, term, x

integer :: n, i

print *, “type value of x and number of terms”

read *,  x, n

sum = x*x/2.0 

term = sum

do i = 1, n

term = term *x*x *(2.0*i-1.0)/(2.0*i+2.0)

sum = sum + term

end do

sum = 1.0 – sum

print *, “value of function is”, sum

end program sum_series

Output of Program

•      x = 0.8       f(x) = 0.6        

type value of x and number of terms

0.8 10

value of function is   0.6000773

type value of x and number of terms

0.8 20

value of function is   0.6000004

type value of x and number of terms

0.8 30

value of function is   0.6000000

•      x = 0.99        f(x)     = 0.14106736..

type value of x and number of terms

0.99 10

value of function is   0.2038259

type value of x and number of terms

0.99 100

value of function is   0.1422002

type value of x and number of terms

0.99 1000

value of function is   0.1410682

Results

•     series converges geometrically

•     possible to find a bound on truncation error

•     number of terms can be selected to get desired bound on error

•     if value of term is too small compared to sum it will be treated as 0.0

•     only a finite number of terms will contribute to sum

–   this number is small if x is not close to 1.0

Integrating Series

•     under certain conditions, possible to integrate a function by integrating each term in its series expansion

•     evaluate integral from –1.0 to 1.0 of sqrt(1.0-x**2)

 

•      integral   = 2 – a0 – a1 – a2 - ……  ,  where

                    a0 = 1/3   and

an = an-1*(2n-1)*(2n+1)/(2n+2)/(2n+3)

 

•      1.5707963268.. is the correct value ( pi/2 )

number of terms

10

value of integral is   1.5755856

number of terms

100

value of integral is   1.5709801

number of terms

1000

value of integral is   1.5708029

Results

•     series does not converge geometrically

•     number of computations is less than in the trapezoidal method

•     increasing number of terms does not reduce error as terms become too small

•     better to add from smallest term to largest

–   more terms can contribute to the sum

Summing Series

•      1.5707963268.. is the correct value ( pi/2 )

number of terms

1000

value of integral is   1.5708020

number of terms

2000

value of integral is   1.5707989

number of terms

4000

value of integral is   1.5707968

Summary

•     complicated mathematical problems can be solved numerically

•     truncation and round-off errors can occur

•     trapezoidal method for evaluating integrals

•     series method for evaluating functions

–   can be used for integration or solving differential equations

–   series summed from smallest to largest