Numerical Computing COL 100 - Introduction to Computer Science - - PowerPoint PPT Presentation

Numerical Computing

COL 100 - Introduction to Computer Science

Department of Computer Science and Engineering Indian Institute of Technology Delhi

Numerical Computing

• Solving mathematical problems numerically

– in terms of NUMBERS, on a computer

• Alternative solution method to difficult

engineering problems

– closed form solution may not be straightforward

• Example problems:

– Integration and differentiation – Solution to equations – Interpolation

Courtesy Prof PR Panda CSE Department IIT Dehi

Numerical Integration

• Computing a definite

integral

• Integral of f(x) between

x=a and x=b is the area under the curve y=f(x) between x=a and b

• This area can be

approximated numerically

y = f(x) x = a x = b

Courtesy Prof PR Panda CSE Department IIT Dehi

Rectangular Rule

• Divide interval [a,b] into n

equal sub-intervals

– length of sub-interval h=(b-a)/n

• In each interval,

approximate f(x) by f(x*)

– f(x*) is the value of f(x) at the mid-point x*of the interval

• Areas of the rectangles are

h.f(x1*), h.f(x2*),...,h.f(xn*) x = a x = b y = f(x) x = a x1* x2* xn* Integral = h.f(x1*) + h.f(x2*) + ... + h.f(xn*) = h. ∑i f(xi*)

Courtesy Prof PR Panda CSE Department IIT Dehi

Program for Numerical Integration

• f is approximated by a

step-function

• This will work no matter

how complex f is

• Accuracy depends on n

– trade-off against computation time

• Only definite integral, no

indefinite integral using this method

def f(x): …. def integrate_midpoint(f,a,b,n): h = (b-a)/n x = a + h/2 sum = 0.0 for i in range(n): sum += f(x) x += h return sum*h

Courtesy Prof PR Panda CSE Department IIT Dehi

Trapezoidal Rule

• Example

https://www.southampton.ac.uk/~fangohr/training/python/pdfs/Python-for-Computational-Science-and-Engineering.pdf

Trapezoidal Rule

• Simple rule

https://www.southampton.ac.uk/~fangohr/training/python/pdfs/Python-for-Computational-Science-and-Engineering.pdf

Trapezoidal Rule

• Example

https://www.southampton.ac.uk/~fangohr/training/python/pdfs/Python-for-Computational-Science-and-Engineering.pdf

Trapezoidal Rule

• Composite Rule

https://www.southampton.ac.uk/~fangohr/training/python/pdfs/Python-for-Computational-Science-and-Engineering.pdf

Trapezoidal Rule

• Divide interval [a,b] into n

equal sub-intervals

– length of sub-interval h=(b-a)/n

• In each interval,

approximate f(x) by line segments with end-points:

– [a,f(a)], [x1,f(x1)],..., [b,f(b)]

• Areas of the trapeziums

are

– 1/2 h (f(a) + f(x1)), – 1/2 h (f(x1) + f(x2)),..., – 1/2 h (f(xn-1) + f(b))

x = b y = f(x) x = a x1 x2 xn-1 Integral = h (f(a)/2 + f(x1) + f(x2) + ... + f(xn-1) + f(b)/2)

Courtesy Prof PR Panda CSE Department IIT Dehi

Trapezoidal Rule

• Divide interval [a,b] into n

equal sub-intervals

– length of sub-interval h=(b-a)/n

• In each interval,

approximate f(x) by line segments with end-points:

– [a,f(a)], [x1,f(x1)],..., [b,f(b)]

• Areas of the trapeziums

are

– 1/2 h (f(a) + f(x1)), – 1/2 h (f(x1) + f(x2)),..., – 1/2 h (f(xn-1) + f(b))

Integral = h (f(a)/2 + f(x1) + f(x2) + ... + f(xn-1) + f(b)/2)

def f(x): …. def integrate_trapezoidal(f,a,b,n): h = (b-a)/n x = a + h sum = (f(a)+f(b))/2 for i in range(1,n): sum += f(x) x += h return sum*h

Courtesy Prof PR Panda CSE Department IIT Dehi

Solving the Equation f(x)=0

• Start with a guess x = x0
• Iteratively refine, attempting to get

closer to the root

• Stop if we detect that we are

– converging: differences between successive x’s get very small – diverging: differences between successive x’s get larger

Courtesy Prof PR Panda CSE Department IIT Dehi

Fixed Point Iteration

• Transform f(x) = 0

algebraically to get

x = g(x) e.g., x2 - 2x +7 = 0 becomes x = 1/2 (x2 + 7)

• Start with arbitrary x = x0
• Compute

– x1 = g(x0), – x2 = g(x1), – ..., – xn = g(xn-1)

• Solution is called fixed point
• f g
• This is a solution to f(x)=0

y = g(x) y = x x = g(x) [solution to f(x)=0]

Courtesy Prof PR Panda CSE Department IIT Dehi

Fixed Point Iteration Example

• Let f(x) = x2 - 3x + 1 = 0
• x = 1/3 (x2 + 1)
• Solutions: x = 2.618, 0.382
• For numerical solution

– let x0 = 1 – x1 = 0.667 – x2 = 0.481 – x3 = 0.411 – x4 = 0.390 – ...converging. Solution!

• Try different starting point

– let x0 = 3 – x1 = 3.333 – x2 = 4.037 – x3 = 5.766 – x4 = 11.415 – ...diverging. Try different start.

y = g(x) y = x 0.382 2.618

Courtesy Prof PR Panda CSE Department IIT Dehi

Bisection Method

• Intermediate value theorem

– if f(a) < 0 and f(b) > 0, then f(x) = 0 somewhere in the interval (a,b)

• Find some a and b for which sign
• f f(x) are different
• Consider c=(a+b)/2
• If f(c) and f(a) are of same sign

– new interval is (c, b) [or (b,c)]

• Else

– new interval is (a,c) [or (c,a)]

• Iterate until interval is small

enough a b c =(a+b)/2 d =(a+c)/2 y=f(x)

Courtesy Prof PR Panda CSE Department IIT Dehi

Newton-Raphson Method

• Another method for solving

f(x)=0

• Assume f ’(x) is continuous
• f ’(x0) = f(x0)/(x0-x1)
• x1 = x0 - f (x0)/f ’(x0)
• x2 = x1 - f (x1)/f ’(x1)
• ...

y = f(x)

x0 x1 x2 x3

xn+1 = xn - f (xn)/f ’(xn)

Courtesy Prof PR Panda CSE Department IIT Dehi

Program for Newton-Raphson Method

• Convergence

detected if new_guess is less than err away from guess

• Fails if tangent slope

is zero

• Signal failure if no

convergence after N iterations

– if failure, try again with different x0

def f(x): def fp(x): def newton(x0, f, fp, error, n): guess = x0 for i in range(n): d = fp(guess) if d == 0: print('There is an error') else: new_guess = guess-f(guess)/d if abs(new_guess - guess) < error: print('successful') return new_guess guess = new_guess print('failed') return 0

Courtesy Prof PR Panda CSE Department IIT Dehi