1 / 33

CS101: Numerical Computing 2

Abhiram Ranade. CS101: Numerical Computing 2. Topics. More examples of for statement Evaluating e x . Hemachandra Numbers Newton Rhapson method for finding roots Defining Functions. e x. We use the series e x = 1 + x + x 2 /2! + x 3 /3! + ...

dior
Télécharger la présentation

CS101: Numerical Computing 2

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Abhiram Ranade CS101: Numerical Computing 2

  2. Topics • More examples of for statement • Evaluating ex. • Hemachandra Numbers • Newton Rhapson method for finding roots • Defining Functions

  3. ex We use the series ex = 1 + x + x2/2! + x3/3! + ... We will evaluate the series to some 10 terms. ith term is xi/i!. How can we calculate it quickly? Direct method: write a loop to calculate xi and also i!, then divide. loops within loop.

  4. Better idea ith term = xi/i! = xi-1/(i-1)! * (x/i) = i-1th term * (x/i) So we calculate ith term from the i-1th term calculated earlier, rather than directly.

  5. Program float x; cin >> x; float term = 1.0; // 0th term of the series. float answer = 0.0; for( int i=1; i <= 10; i++){ answer += term; term = term * x / i;} cout << answer;

  6. Notes • You must put the following comment to your program: /* After ith iteration, i terms of the series have been summed up */ • Check that this is true! - by inspection, - Add a statement to print after each iteration and check.

  7. Program 2 float x; cin >> x; float term = 1; // 0th term of the series. float answer = 0; for( int i=1; term > 0.0001; i++){ answer += term; term = term * x / i;} cout << answer;

  8. Remarks • If we calculate xi/i! independently in each iteration, we will need about 2i arithmetic operations in each iteration. • In our program, we just need 2 in each iteration. • Idea of reusing what was calculated in the previous iteration is very important, and also very commonly useful.

  9. Hemachandra’s Problem (12th century AD) • Suppose I have to build a wall of length 8 feet. I have bricks 2 feet long and also 1 foot long. In how many ways I can lay the bricks so that I fill the 8 feet? • Possiblities: 2,2,2,2; 1,1,1,1,1,1,1,1; 2,2,2,1,1 ....

  10. Hemachandra’s Actual Problem • Suppose I am designing a poetic meter with 8 beats. The meter is made of short syllables and long syllables. Short syllable = 1 beat, Long syllable = 2 beats. How many ways are there of filling 8 beats? • Example of a poetic meter Aa ji chya Ja wa lii gha dya l ka sa lay aa he ch mat kaa ri k ya kun den du tu sha r ha r dha va la ya shubh r vas tra vru ta l l l s s l s l s s s l l l s l l s s

  11. Hemachandra’s Solution • “By the method of Pingala, it is enough to observe that the last beat is long or short” • Pingala: mathematician/poet from 500 A.D. • Hemachandra is giving credit to someone who lived 1500 years before him!! • Copy but give credit..

  12. Hemachandra’s solution contd. S : Class of 8 beat patterns with short last beat. L : Class of 8 beat patterns with long last beat. Each 8 beat pattern is in class L or class S S = all 7 beat patterns + long beat appended. | class S | = Number of patterns with 7 beats | class L | = Number of patterns with 6 beats |8 beat patterns| = |class S| + |class L| = |7 beat patterns| + |6 beat patterns|

  13. Algebraically.. Hn = number of patterns with n beats H8 = H7 + H6 In general Hn = Hn-1 + Hn-2 Does this help us to compute H8? We need to know H7, H6, for which we need H5, ...

  14. Algorithm Idea • H1 = number of patterns with 1 beat = 1 • H2 = Number with 2 beats = 2 ...{SS, L} • H3 = H2 + H1 = 2 + 1 = 3 • H4 = H3 + H2 = 3 + 2 = 5 • H5 = H4 + H3 = 5 + 3 = 8 • H6 = H5 + H4 = 8 + 5 = 13 • H7 = H6 + H5 = 13 + 8 = 21 • H8 = H7 + H6 = 21 + 13 = 34 ...

  15. Program to compute Hn int n; cin >> n; // which number to compute int hprev = 1, hcurrent = 2; for(int i=3; i <= n; i++){ hnext = hprev + hcurrent; hprev = hcurrent; // prepare for next iteration hcurrent = hnext; } cout << hnext;

  16. Code is tricky! • Need a comment: /* At the begining of an iteration hcurrent = Hi-1 write h[i] if you like. hprev = Hi-2 where i is the value of variable i */ • Can you prove this? • Will mathematical induction help? • Proving this is enough -- hnext = hprev + hcurrent -- hence correct answer will be generated.

  17. Proof by induction • Base case: At the beginning of the first iteration is this true? Yes, i will have value 3, and hprev = 1 = H1, hcurrent = 2 = H2 • Suppose it is true at some later iteration, when i has value v >= 3. By induction hypothesis hprev, hcurrent have values Hv-1,Hv-2. • The first statement hnext = hprev + hcurrent makes hnext = Hv-1 + Hv-2 = Hv. After this the statement hprev = hcurrent makes hprev = Hv-1. The next statement hcurrent = hnext makes hcurrent=Hv. In the next iteration i will have value v+1. But hprev,hcurrent will have exactly the right values!

  18. On Hemachandra Numbers • Series is very interesting. • Number of petals in many flowers. • Ratio of consecutive terms tends to a limit. • What are these numbers more commonly known as? • Hemachandra lived before Fibonacci. • Mathematics from poetry!

  19. Newton Raphson method Method to find the root of f(x), i.e. x s.t. f(x)=0. • Method works if: • f(x) and f '(x) can be easily calculated. • A good initial guess is available. • Example: To find square root of k. • use f(x) = x2 - k. f ‘ (x) = 2x. • f(x), f ‘ (x) can be calculated easily. 2,3 arithmetic ops. • Initial guess x0 = 1 always works! can be proved.

  20. How to get better xi+1 given xi Point A =(xi,0) known. Calculate f(xi ). Point B=(xi,f(xi)) known B Approximate f by tangent C= intercept on x axis C=(xi+1,0) C A xi+1 xi f(x) f ‘ (xi) = AB/AC = f(xi)/(xi - xi+1)  xi+1 = (xi- f(xi)/f ‘ (xi))

  21. Square root of k xi+1 = (xi- f(xi)/f ‘ (xi)) f(x) = x2 - k, f ‘ (x) = 2x xi+1 = xi - (xi2 - k)/(2xi) = (xi + k/xi)/2 Starting with x0=1, we compute x1, then x2, then... We can get as close to sqrt(k) as required. Proof not part of the course.

  22. Code float k; cin >> k; float xi=1; // Initial guess. Known to work. for(int i=0; i < 10; i++){ // 10 iterations xi = (xi + k/xi)/2; } cout << xi;

  23. Another way float xi, k; cin >> k; for( xi = 1 ; // Initial guess. Known to work. xi*xi – k > 0.001 || k - xi*xi > 0.001 ; // until error in the square is at most 0.001 xi = (xi + k/xi)/2); cout << xi;

  24. Yet Another way float k; cin >> k; float xi=1; while(xi*xi – k > 0.001 || k - xi*xi > 0.001){ xi = (xi + k/xi)/2 ; } cout << xi; }

  25. While statement • while (condition) { loop body} • check condition, then execute loop body if true. Repeat. • If loop body is a single statement, then need not use { }. Always putting braces is recommended; if you insert a statement, you may forget to put them, so do it at the beginning. • True for other statements also: for/repeat/if.

  26. For vs. while • If there is a “control” variable with initial value, update rule, and whose value distinctly defines each loop iteration, use for. • If loop executes fixed number of times, use for. • Personal taste. • “break” and “continue” : read from your book.

  27. Procedures which return values • Polygon(sides, length) draws. • But it will be nice to have a procedure sqrt which can be called as y = sqrt(2.56); // should return 1.6 which would suspend the current program, compute sqrt(k), send back the value, and resume. • This is possible!

  28. Procedure for sqrt float sqrt(float k){ float xi=1; while(xi*xi – k > 0.001 || k - xi*xi > 0.001){ xi = (xi + k/xi)/2 ; } return xi; }

  29. Notes • First word tells the type of the value being returned. • Note the return statement: this says what value is to be sent back. • Other nomenclature is same. Parameters (k), argument (2.56 in our example). • “y = sqrt(2.56);” is a “call”.

  30. Execution Mechanism • Almost identical to what was described earlier. • Calling program suspends. • Area constructed for sqrt procedure. Argument value (2.56) copied to k from main program. • Sqrt program executes in its own area. • After it finishes, result (1.6 if it works) is copied back to calling program. The result replaces the call. Area in which sqrt executed is destroyed. • Calling program resumes.

  31. Notes • In C++ procedures are called functions. • Uniform syntax for functions that do not return values. So far we have used: “ procedure Polygon(...)...” Instead in C++ it is more common to write “ void Polygon(...) ...” • “Procedure” was a keyword defined only for this course, stop using it from now on.

  32. Math library You didnt actually need to write the sqrt function. Just add the line: #include <cmath> at the top of your file. This allows you to use in your program the mathematical functions in the cmath library. sqrt is in the library. So are functiones such as sine, cos, tan, ...

  33. Homework • Write a program to calculate the cube root using Newton Raphson method. • Check how many iterations are needed to get good answers. Should be very few. • Use sqrt function to draw a right angled isoceles triangle.

More Related