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 y = f(x) • Integral of f(x) between x=a and x=b is the area under the curve y=f(x) between x=a and b x = a x = b • This area can be approximated numerically Courtesy Prof PR Panda CSE Department IIT Dehi

Rectangular Rule y = f(x) • 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 x = a x = b x = a mid-point x*of the interval • Areas of the rectangles are h.f(x 1 *), h.f(x 2 *),...,h.f(x n *) x 1 * x 2 * x n * Integral = h.f(x 1 *) + h.f(x 2 *) + ... + h.f(x n *) = h. ∑ i f(x i *) Courtesy Prof PR Panda CSE Department IIT Dehi

Program for Numerical Integration • f is approximated by a step-function def f(x): …. • This will work no matter def integrate_midpoint(f,a,b,n): how complex f is h = (b-a)/n • Accuracy depends on n x = a + h/2 sum = 0.0 – trade-off against for i in range(n): computation time sum += f(x) • Only definite integral, no x += h indefinite integral using return sum*h this method 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 y = f(x) – length of sub-interval h=(b-a)/n • In each interval, approximate f(x) by line segments with end-points: – [a,f(a)], [x 1 ,f(x 1 )],..., [b,f(b)] • Areas of the trapeziums x = a x = b are – 1/2 h (f(a) + f(x 1 )), – 1/2 h (f(x 1 ) + f(x 2 )),..., x 1 x 2 x n-1 – 1/2 h (f(x n-1 ) + f(b)) Integral = h (f(a)/2 + f(x 1 ) + f(x 2 ) + ... + f(x n-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 def f(x): …. h=(b-a)/n • In each interval, def integrate_trapezoidal(f,a,b,n): approximate f(x) by line h = (b-a)/n segments with end-points: x = a + h – [a,f(a)], [x 1 ,f(x 1 )],..., [b,f(b)] sum = (f(a)+f(b))/2 • Areas of the trapeziums for i in range(1,n): are sum += f(x) x += h – 1/2 h (f(a) + f(x 1 )), return sum*h – 1/2 h (f(x 1 ) + f(x 2 )),..., – 1/2 h (f(x n-1 ) + f(b)) Integral = h (f(a)/2 + f(x 1 ) + f(x 2 ) + ... + f(x n-1 ) + f(b)/2) Courtesy Prof PR Panda CSE Department IIT Dehi

Solving the Equation f(x)=0 • Start with a guess x = x 0 • 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 y = g(x) x = g(x) e.g., x 2 - 2x +7 = 0 becomes x = 1/2 (x 2 + 7) • Start with arbitrary x = x 0 • Compute – x 1 = g(x 0 ), – x 2 = g(x 1 ), y = x – ..., – x n = g(x n-1 ) • Solution is called fixed point x = g(x) of g [solution to f(x)=0] • This is a solution to f(x)=0 Courtesy Prof PR Panda CSE Department IIT Dehi

Fixed Point Iteration Example • Let f(x) = x 2 - 3x + 1 = 0 • x = 1/3 (x 2 + 1) • Solutions: x = 2.618, 0.382 • For numerical solution – let x 0 = 1 y = g(x) – x1 = 0.667 y = x – x2 = 0.481 – x3 = 0.411 – x4 = 0.390 – ...converging. Solution! • Try different starting point – let x 0 = 3 – x1 = 3.333 – x2 = 4.037 – x3 = 5.766 – x4 = 11.415 – ...diverging. Try different start. 2.618 0.382 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 of f(x) are different y=f(x) • Consider c=(a+b)/2 • If f(c) and f(a) are of same sign – new interval is (c, b) [or (b,c)] • Else a – new interval is (a,c) [or (c,a)] c • Iterate until interval is small b =(a+b)/2 d enough =(a+c)/2 Courtesy Prof PR Panda CSE Department IIT Dehi

Newton-Raphson Method • Another method for solving f(x)=0 y = f(x) • Assume f ’(x) is continuous • f ’(x 0 ) = f(x 0 )/(x 0 -x 1 ) • x 1 = x 0 - f (x 0 )/f ’(x 0 ) • x 2 = x 1 - f (x 1 )/f ’(x 1 ) • ... x n+1 = x n - f (x n )/f ’(x n ) x 3 x 2 x 1 x 0 Courtesy Prof PR Panda CSE Department IIT Dehi

Program for Newton-Raphson Method • Convergence def f(x): detected if def fp(x): new_guess is less than err away from def newton(x0, f, fp, error, n): guess = x0 guess for i in range(n): • Fails if tangent slope d = fp(guess) if d == 0: print('There is an error') is zero else: new_guess = guess-f(guess)/d • Signal failure if no if abs(new_guess - guess) < error: convergence after N print('successful') iterations return new_guess guess = new_guess – if failure, try again print('failed') with different x0 return 0 Courtesy Prof PR Panda CSE Department IIT Dehi