Numerical Methods

1. Numerical Methods Applications of Loops: The power of MATLAB Mathematics + Coding

2. Numerical Methods Numerical methods are techniques used to find quantitative solutions to mathematical problems. Many numerical methods are iterative solutions– this means that the technique is applied over and over, gradually getting closer (i.e. converging) on the final answer. Iterative solutions are complete when the answer from the current iteration is not significantly different from the answer of the previous iteration. “Significantly different” is determined by the program or user – for example, it might mean that there is less than 1% change or it might mean less than 0.00001% change.

3. Root Finding • If we can find where f(x) = 0, then we can find where f(x) equals any value. • To find where f(x) = -2, find the roots of f(x)+2 • Essentially, we can shift any function up or down and solve for the roots, so all value finding problems can essentially be adjusted to be root finding problems. • We can turn any equation into a root finding problem.

4. Bracketing Methods • There are 2 basic “bracketing methods” (only 1 of which we’ll discuss). • They both rely on the same basic concept • If you know that f(x1) > 0 and f(x2) < 0, then a root must occur between x1 and x2.

5. More ExamplesES206 Fluid Dynamics Solve for f The Colebrook Equation

6. ES206 – Fluid Dynamics “WOW – that’s hard to solve!” How about this instead: “By hand, find the square root of 127.3” Both problems use the same technique to solve: the Bisection Method.

7. Bisection Method • For the Colebrook Equation, we know that this equation must hold: • Re-arrange and you have a roots problem. Guess for f and if the equation is not satisfied, adjust f and try again. • This is the same as solving for 127.3 = X2. We made it X2 -127.3 = 0 and found the roots.

8. Bisection Method • The Bisection Method is a very common technique for solving mathematical problems numerically. • Most people have done this even if they didn’t know the name. • The technique involves taking a guess for the desired variable, then applying that value to see how close the result is to a known value. • Based upon the result of the previous step, a better “guess” is made and the test repeated.

9. Bisection Method Typical human method: Find the square root of 127.3 • First guess: 11 (since we know that 112 is 121) • 11*11 = 121  too low (looking for 127.3) • Adjust: 12 is the new guess • 12*12 = 144  too high • Adjust: 11.5 (“bisect” the difference) • 11.5*11.5 = 132.25  too high • Adjust: 11.25 • 11.25 * 11.25 = 126.5625  too low • etc.

10. Bisection Method Let’s use the programmer’s method – an algorithm: • Define x1 so that f(x1) is greater than 0. • Define x2 so that f(x2) is less than 0. • Evaluate f((x1+x2)/2) • If it is positive, set x1 to (x1+x2)/2 • If is it negative, set x2 to (x1+x2)/2 • Repeat until you’re happy with x1-x2

11. Pros/Cons with Bisection • Cons • Requires f(x1) < 0 and f(x2) > 0. • CAN be hard to find • Can require A LOT of iterations. • Pros • GUARANTEED to work

12. Newton’s Method X is the root we’re trying to find – the value where the function becomes 0 Xn is the initial guess (hopefully, a reasonably-close guess!) Xn+1 is the next approximation using the tangent line Credit to: http://en.wikipedia.org/wiki/Newton’s_method

13. The Big Picture (1/4) 1st: guess! xn

14. The Big Picture (1/4) 2nd: Go hit the actual curve At this (xn, yn), calculate the slope of the curve using the derivatives: m = f’(xn)

15. The Big Picture (1/4) 3rd: Using the slope m, the coordinates (Xn , Yn), calculate Xn+1 : Remember: m was obtained from the derivative – only Xn+1 is unknown. What is Yn+1?

16. The Big Picture (2/4) TOO FAR Guess again!

17. The Big Picture (2/4) Iterate!

18. The Big Picture (2/4) TOO FAR Guess again!

19. The Big Picture (3/4) – Eventually… CLOSE ENOUGH! New estimate Old estimate Stop the loop, when there is not much difference between the guess and the new value!

20. Newton’s Method Prepare a MATLAB program that uses Newton’s Method to find the roots of the equation: y(x) = x3 + 2x2 – 300/x

21. Newton’s Method f(x) = x3 + 2x2 – 300/x What does it seem the answer is?

22. Newton’s Method Closer look… Actual Root!

23. Ex 1: algorithm, skeleton % Provide an initial guess for the x-intercept % Repeatuntil close enough: compare most recent x-intercept to previous x-intercept % Find the slopeof the tangent line at the point (requires derivative) % Using the slope formula, find the x intercept of the tangent line

24. Newton’s Method % Define initial difference as very large (force loop to start) % Repeat as long as change in x is large enough % Make the “new x” act as the current x % Find slope of tangent line using f’(x) % Compute y-value % Use tangent slope to find the next value: % two points on line are the current point (x, y) % and the x-intercept (x_new, 0)

25. % Define initial difference as very large • % (force loop to start) • x_n = 1000; • x_np1 = 5; % initial guess • % Repeat as long as change in x is large enough • while (abs(x_np1 – x_n) > 0.0000001) • % x_np1 now becomes the new guess • x_n = x_np1; • % Find slope of tangent line f'(x) at x_n • m = 3*x_n^2 + 4*x_n + 300/(x_n^2); • % Compute y-value • y = x_n^3 + 2*x_n^2 - 300/x_n; • % Use tangent slope to find the next value: • % two points on line: the current point (x_n, y_n) • % and the x-intercept (x_np1, 0) • x_np1 = ((0 - y) / m) + x_n; • end Xn+1 Xn+1 TOO FAR ? Xn Xn

26. Newton’s Method Add fprintf’s and see the results: x y ------------------- ------------------- 5.000000000000000 115.000000000000000 3.925233644859813 14.864224007996924 3.742613907949879 0.279824373263295 3.739045154311932 0.000095471294827 3.739043935883257 0.000000000011084 Root of f(x) = x3+2x2–300/x is approximately 3.7390439

27. Real-World • Numerical Methods are used in CFD, where convergence can take weeks/months to complete (get “close enough”)... • while loops are generally used, since the number of iteration depends on how close the result should be, i.e. the precision. • Practice! Code this, and usea scientific calculator to see how close the results are.

28. Cons • Cons • Divide by 0 Error • What happens if x is a maximum or minimum – the slope of the tangent line is 0. • Divergence • What if our tangent line leads us away from our root. • Cycling • f(x)=x3-x-3 • f’(x)=3x2-1 • x0=0

29. Numerical Methods Secant Method: Find the Derivative

30. Derivatives, intro. Another numerical method involves finding the value of the derivative at a point on a curve. Basically, this method uses secant lines as approximations to the tangent line. This method is used when the actual symbolic derivative is hard to find, the formula is simply not needed, or the function is prone to change as the main program does its job.

31. Derivatives, big picture What is the slope at this point?

32. Derivatives (1/5) Choose a and b to be equidistant from that point, calculate slope 1 Slope 1 range X range/2 range/2

33. Derivatives (2/5) Reduce the range for the second iteration. Slope 1 x loX hiX Zoom in on x

34. Derivatives (3/5) Slope 1 x loX hiX Find new secant line.

35. Derivatives (3/5) Slope 1 Slope 2 hiY loY x loX hiX Calculate new slope

36. Derivatives (3/5) Slope 1 Slope 2 Is slope 2 “significantly different” from slope 1? x loX hiX Assume they are different… let’s iterate!...

37. Derivatives (4/5) Slope 1 Slope 2 Slope3 x loX hiX Zoom in again..

38. Derivatives (4/5) Slope 1 Slope 2 Slope3 Is Slope2 very different from Slope3? Iterate while yes! x loX hiX Zoom in again..

39. (5/5) Eventually… Stop the loop when the slopes are no longer “significantly different” Remember: “significantly different” means • maybe a 1% difference? • maybe a 0.1% difference? • maybe a 0.000001% difference?

40. Derivatives - Algorithm % Define initial difference as very large % (force loop to start) % Specify starting range % Repeat as long as slope has changed “a lot” % Cut range in half % Compute the new low x value % Compute the new high x value % Compute the new low y value % Compute the new high y value % Compute the slope of the secant line % end of loop Let’s do this for f(x) = sin(x) + 3/sqrt(x)

41. Code % Define initial difference as very large % (force loop to start) m = 1; oldm = 2; % Specify starting range and desired point range = 10; x=5; ctr=0; % Repeat as long as slope has changed a lot while (ctr<3 || (abs((m-oldm)/oldm) > 0.00000001)) % Cut range in half range = range / 2; % Compute the new low x value lox = x – range ; % Compute the new high x value hix = x + range ; % Compute the new low y value loy = sin(lox) + 3/sqrt(lox); % Compute the new high y value hiy = sin(hix) + 3/sqrt(hix); % Save the old slope oldm = m; % Compute the slope of the secant line m = (hiy – loy) / (hix – lox); % increment counter ctr = ctr + 1; end% end of loop (or have just started)

42. Derivatives - the end Add some fprintf’s and… LOW x HIGH x Slope Difference in slope 1: 0.000000000000000, 10.000000000000000 -Inf -Inf 2: 2.500000000000000, 7.500000000000000 -0.092478729683983 Inf 3: 3.750000000000000, 6.250000000000000 0.075675505484728 0.168154235168711 4: 4.375000000000000, 5.625000000000000 0.130061340134248 0.054385834649520 5: 4.687500000000000, 5.312500000000000 0.144575140383541 0.014513800249293 6: 4.843750000000000, 5.156250000000000 0.148263340290110 0.003688199906569 7: 4.921875000000000, 5.078125000000000 0.149189163013407 0.000925822723297 8: 4.960937500000000, 5.039062500000000 0.149420855093148 0.000231692079740 9: 4.980468750000000, 5.019531250000000 0.149478792897432 0.000057937804284 10: 4.990234375000000, 5.009765625000000 0.149493278272689 0.000014485375257 11: 4.995117187500000, 5.004882812500000 0.149496899674250 0.000003621401561 12: 4.997558593750000, 5.002441406250000 0.149497805028204 0.000000905353954 13: 4.998779296875000, 5.001220703125000 0.149498031366966 0.000000226338761 14: 4.999389648437500, 5.000610351562500 0.149498087951815 0.000000056584850 15: 4.999694824218750, 5.000305175781250 0.149498102098005 0.000000014146190 16: 4.999847412109375, 5.000152587890625 0.149498105634484 0.000000003536479 17: 4.999923706054688, 5.000076293945313 0.149498106518877 0.000000000884393 Slope of f(x) = sin(x) + 3/sqrt(x) at x=5is approximately 0.1494981

43. Integration

44. Integration • Some functions can’t be integrated symbolically • Other functions are just REALLY ugly to do it

45. Trapezoidal Rule • Approximate the curve with an easy to integrate first order polynomial - This ends up being a trapezoid.

47. Trapezoidal Rule • Trapezoidal Rule works out as:

48. Composite Trapezoidal • Can we just combine multiple trapezoids? • Yes - Instead of 1 large trapezoid, we’ll use multiple smaller trapezoids and add up their areas (and their errors)