CS 101 Numerical Computing 2 Abhiram Ranade Topics

CS 101: Numerical Computing 2 Abhiram Ranade

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

ex We use the series ex = 1 + x 2/2! + x 3/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.

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

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

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.

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

Remarks If we calculate xi/i! independently in each iteration, we will need about 2 i 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.

Hemachandra’s Problem (12 th 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; 2, 2, 2, 1, 1. .

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 ya kun den du tu sha r l l l s s l s ch mat kaa ri k ha r dha va la ya shubh r vas tra vru ta l s s s l l l s s

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. .

Algebraically. . Hn = number of patterns with n beats H 8 = H 7 + H 6 In general Hn = Hn-1 + Hn-2 Does this help us to compute H 8? We need to know H 7, H 6, for which we need H 5, . . .

Algorithm Idea H 1 = number of patterns with 1 beat = 1 H 2 = Number with 2 beats = 2. . . {SS, L} H 3 = H 2 + H 1 = 2 + 1 = 3 H 4 = H 3 + H 2 = 3 + 2 = 5 H 5 = H 4 + H 3 = 5 + 3 = 8 H 6 = H 5 + H 4 = 8 + 5 = 13 H 7 = H 6 + H 5 = 13 + 8 = 21 H 8 = H 7 + H 6 = 21 + 13 = 34. . .

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;

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.

Proof by induction Base case: At the beginning of the first iteration is this true? Yes, i will have value 3, and hprev = 1 = H 1, hcurrent = 2 = H 2 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!

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!

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) = x 2 - k. f ‘ (x) = 2 x. f(x), f ‘ (x) can be calculated easily. 2, 3 arithmetic ops. Initial guess x 0 = 1 always works! can be proved.

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

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

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

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;

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; }

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.

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.

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!

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; }

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”.

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.

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.

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, . . .

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.