CS3911 Intro to Numerical Methods with Fortran Final Solutions 1 CS3911 Introduction to Numerical Methods with Fortran by yyy55749

VIEWS: 104 PAGES: 10

• pg 1
```									CS3911 Intro. to Numerical Methods with Fortran Final Solutions                                  1

CS3911 Introduction to Numerical Methods with
Fortran
Final Exam Solutions
1. Accuracy and Reliability

(a) [20 points] Let be a very small positive number (i.e., ≈ 0). Do the following two
problems: (1) solve the following system of linear equations with and without partial
pivoting, assuming you use ﬁnite precision arithmetics; and (2) ﬁnd the major reason
or reasons that can explain the diﬀerence(s) between the two solutions. You have
to clearly state your ﬁndings with a convincing argument. Just stating a
“reason” such as “it is because of cancelation” or “overﬂow” will receive zero
point.
                                   
1           x               1
                                   
           ·         =             
1 1             y               2

Answer: Without pivoting, one multiplies −1/ to the ﬁrst equation and adds the result
to the second. This yields the following:
                                                        
1               x                   1
                                                        
                       ·         =                      
0 1 − 1/                    y               2 − 1/

Since is small, 1/ is large. As a result, -1/ is the dominating term in both 1 − 1/
and 2 − 1/ . In other words, 1 − 1/ ≈ −1/ and 2 − 1/ ≈ −1/ . Hence, the equation
(1 − 1/ )y = 2 − 1/ numerically becomes
1           1
−       y≈−

Therefore, we have y = 1. Plugging y = 1 into the ﬁrst equation x + y = 1 yields x = 0.
This is certainly incorrect because x = 0 and y = 1 do not satisfy the second equation
x + y = 2.
With pivoting, swapping the equations transforms the system to the following:
                                   
1 1             x               2
                                   
           ·         =             
1           y               1

Multiplying the ﬁrst equation by − and adding the result to the second yields:
                                                
1       1               x                   2
                                                
                   ·         =                  
0 1−                    y               1−2

Since is small, it does not contribute much to 1 − and 1 − 2 . As a result, 1 − ≈ 1
and 1 − 2 ≈ 1, and the second equation (1 − )y = 1 − 2 numerically becomes y = 1.
Plugging y = 1 into the ﬁrst equation gives x = 1.
Obviously, x = y = 1 is the correct numerical solution if is very small. The reason for
the non-pivoting solution to go wrong is rounding error in computing 1 − 1/ and 2 − 1/
if 1/ is very large. In this case, rounding error makes both terms nearly equal to −1/ .
CS3911 Intro. to Numerical Methods with Fortran Final Solutions                                    2

Consequently, we have y = 1 which is still correct. But, in solving x from x + y = 1
with y = 1, we have x = 0/ = 0.
It does not require an extremely small to cause this problem. For example, on a
7-signiﬁcant-digit computer, if = 0.00000001 we have 1/ = 100, 000, 000. Then,
1 − 1/ = −99, 999, 999 and 2 − 1/ = −99, 999, 998. Both would be rounded to 7-digit
and the result is 0.1 × 109.
This example clearly illustrates that pivoting is necessary.

2. Systems of Linear Equations

(a) [10 points] Answer the following two questions: (1) deﬁne the meaning of a diagonal
dominant matrix A = [ai,j ]n×n ; and (2) what is the advantage(s) for a matrix to be
diagonal dominant?
Answer: A n × n matrix A = [ai,j ] is diagonal dominant if for every row i the following
condition holds:
n
|ai,j | < |ai,i |
j=1,j=i

In other words, for each row, the sum of the absolute values of all oﬀ-diagonal entries is
less than the absolute value of the diagonal element.
If a matrix is diagonal dominant, Gauss elimination, LU-decomposition included, does
not need pivoting, and Jacobi’s and Gauss-Seidel’s methods always converge.

3. Interpolation and Approximation

(a) [10 points] Given a polynomial of degree n as follows,

Pn (x) = a0 + a1 x + a2x2 + a3x3 + · · · + an xn

suggest an eﬃcient way to evaluate Pn (x) with n multiplications. You should provide
an algorithm and discuss its correctness and complexity. A method that does
not achieve n multiplications receives zero point.
Answer: The given polynomial can be rewritten in a nested form as follows:

Pn (x) = a0 + (a1 + (a2 + (a3 + (· · · (an−1 + an x)x) · · ·)x)x)

For example, P4 (x) = a0 + a1 x + a2 x2 + a3 x3 + a4 x4 = a0 + (a1 + (a2 + (a3 + a4x)x)x)x.
In this way, only n multiplications are needed (i.e., O(n)). On the other hand, since
ai xi requires i multiplications (i.e., i − 1 for xi and 1 for ai (xi)), the total number of
multiplications with a direct evaluation of the power form is 0+1+2+· · ·+n = n(n+1)/2
(i.e., O(n2)).
The following is a possible algorithm, which assumes that coeﬃcients ai ’s are in array
a( ), x is in variable x, and Px has the result:
Px = a(n)
DO i = n-1, 0, -1
Px = a(i) + Px*x
END DO
Since this DO loop iterates n times, each of which uses exactly one multiplication, the
order of complexity is O(n).
One may suggest the following O(n) implementation:
CS3911 Intro. to Numerical Methods with Fortran Final Solutions                                       3
Px    = a(0)
Power = x
DO i = 1, n
Px     = Px + a(i)*Power
Power = Power*x
END DO
Although it is an O(n) method, it uses two multiplications per iteration and the total
number of multiplications is 2n. This is certainly slower than the nested form.

(b) [20 points] Suppose ﬁve data points (x0 , y0) = (−2, 1), (x1, y1 ) = (−1, 1), (x2 , y2) =
(0, 1), (x3, y3) = (1, 1) and (x4, y4) = (2, 25) are given. Find the Newton interpolating
polynomial for the given data points with the divided diﬀerence method. You should
show all computation steps. Only providing an answer and/or using a wrong
Answer: The Newton interpolating polynomial for ﬁve data points is

P4 (x) = f [x0] + f [x0 , x1](x − x0) + f [x0 , x1, x2](x − x0)(x − x1 )
+ f [x0, x1, x2, x3](x − x0)(x − x1 )(x − x2 )
+ f [x0, x1, x2, x3, x4](x − x0)(x − x1 )(x − x2)(x − x3 )

The divided diﬀerence table is shown below:

xi   f[xi ]
−2       1

−1        1    →      0

0       1    →      0    →      0

1       1    →      0    →      0    →     0

2      25    →     24    →    12     →     4    →    1

Therefore, we have f [x0 ] = 1, f [x0, x1] = 0, f [x0 , x1, x2] = 0, f [x0, x1, x2, x3] = 0 and
f [x0, x1, x2, x3, x4] = 1, and the interpolating polynomial is P4 (x) = 1 + (x + 2)(x +
1)x(x − 1).

4. Numerical Diﬀerentiation:

(a) [10 points] Use the 3-point backward diﬀerence method to compute the derivative of
f (x) = ex sin(x) at x = 1 with ∆ = 0.01. You should show all computation steps.
Only providing an answer and/or using a wrong method receives zero point.
Answer: The 3-point backward diﬀerence method uses the following formula:
1
f (x) ≈      (3fi − 4fi−1 + fi−2 )
2∆
Since ∆ = 0.01, fi = f (1) = 2.287355, fi−1 = f (1 − 0.01) = 2.249942 and fi−2 =
f (1 − 0.02) = 2.212824, we have
1
f (1) ≈            (3 × 2.287355 − 4 × 2.249942 + 2.212824) = 3.756102
2 × 0.01
CS3911 Intro. to Numerical Methods with Fortran Final Solutions                                         4

(b) [20 points] Use Richardson’s extrapolation method to compute the derivative of f (x) =
ex sin(x) at x = 1 with initial ∆ = 0.1. You should carry out the extrapolation steps until
the computation for ∆ = 0.025 completes, and clearly indicate the desired answer.
You should show all computation steps. Only providing an answer and/or
using a wrong method receives zero point.
Answer: Column 0 (i.e., d0,0, d1,0, d2,0, etc) uses 2-point central diﬀerence. Since the
starting ∆ is 0.1 (i.e., row 0), the ∆’s for row 1 and row 2 are 0.05 = 0.1/2 and 0.025 =
0.5/2, respectively. The following shows d0,0, d1,0 and d2,0 using central diﬀerence with
∆ being 0.1, 0.05 and 0.025:

f (1 + 0.1) − f (1 − 0.1)
d0,0 =                                   = 3.753307857
2 × 0.1
f (1 + 0.05) − f (1 − 0.05)
d1,0 =                                     = 3.755366227
2 × 0.05
f (1 + 0.025) − f (1 − 0.025)
d2,0 =                                       = 3.75587862
2 × 0.025
Richardson’s extrapolation uses the following formula:
di,j − di−1,j
di,j+1 = di,j +
4j+1 − 1
The following computes the desired result in d2,2:

Row    ∆             Column 0                           Column 1                 Column 2
0     0.1      d0,0 = 3.753307857

1     0.05     d1,0 = 3.755366227        →     d1,1 = 3.7560523567

2     0.025    d2,0 = 3.75587862         →     d2.1 = 3.7560494227    →     d2,2 = 3.756049227

Therefore, the extrapolated derivative approximation is d2,2 = 3.756049227.

5. Numerical Integration: The following integral has an exact solution:
1
ex dx = e1 − e0 = e − 1 = 1.71828 . . .
0

(a) [20 points] Use the iterative 3-point Simpson method to compute the above integration.
You should start with ∆ = 0.5, update the result in each iteration, and stop after the
result of ∆ = 0.125 is obtained. You should show all computation steps. Only
Answer: The 3-point Simpson method uses the following formula:
m             m
∆
(f0 + f2m ) + 4     f2i−1 + 2     f2i
3                 i=1           i=1

The table below shows the computed results. In the table, “Previous Sum” is the sum
of terms computed in the previous iteration, which becomes the sum of “even” terms in
this iteration (i.e., m f2i ), and “This Odd Sum” is the sum of newly introduced terms
i=1
in this iteration (i.e., m f2i−1 ). The sum of “Previous Sum” and “This Odd Sum”
i=1
CS3911 Intro. to Numerical Methods with Fortran Final Solutions                                        5

becomes the “Previous Sum” for the next iteration. The desired result is 1.718284155.

∆              0.5                0.25          0.125
f (0) + f (1)     3.718281828        3.718281828    3.718281828
0                      1
0.125                                               1.133148453
0.25                                 1.284025417
0.375                                               1.454991415
0.5               1.648721271
0.625                                               1.868245957
0.75                                 2.117000017
0.825                                               2.398875294
1                 2.718281828
Previous Sum           0             1.648721271    5.649746705
This Odd Sum      1.648721271        3.401025434    6.855261119
Integration       1.718861152        1.718318842    1.718284155

(b) [20 points] Use Romberg’s method to compute the above integration and use the trape-
zoid method for column 0. You should carry out all steps until the computation for
∆ = 0.25 completes, and clearly indicate the desired answer. You should show all
computation steps. Only providing an answer and/or using a wrong method
Answer: Column 0 (i.e., r0,0, r1,0, r2,0, etc) uses the trapezoid method with 1, 2 and 4
subintervals. Since the starting interval is [0, 1], the length of subintervals for rows 0, 1
and 2 are 1, 0.5, and 0.25, respectively. The following shows r0,0, r1,0 and r2,0 using the
trapezoid method with ∆ = 0.1, 0.05 and 0.025:

f (0) + f (1)
r0,0 = ∆                   = 1.859149914
2
f (0) + f (1)
r1,0    = ∆               + f (0.5) = 1.753931092
2
f (0) + f (1)
r2,0    = ∆               + f (0.25) + f (0.5) + f (0.75) = 1.727221905
2
Romberg’s method uses the formula identical to that of Richardson’s extrapolation:
ri,j − ri−1,j
ri,j+1 = ri,j +
4j+1 − 1
The following computes the desired result in r2,2:

Row     ∆             Column 0                     Column 1                Column 2
0      1        r0,0 = 1.859149914

1     0.5      r1,0 = 1.753931092   →     r1,1 = 1.718861150

2     0.25     r2,0 = 1.727221905   →     r2.1 = 1.718318843    →    r2,2 = 1.718282689

Therefore, Romberg’s method yields r2,2 = 1.718282689.
CS3911 Intro. to Numerical Methods with Fortran Final Solutions                                 6

6. ODE Initial Value Problems Suppose you are given an ODE y = xy with initial value
y(0) = 1.

(a) [10 points] Use Euler’s method to solve the given equation on [0, 4] with ∆ = 1, and
ﬁll your results in the following table.
Answer: Euler’s method uses y(x + ∆) = y(x) + ∆f (x, y) to compute the next y. The
following table shows the computed results:

Iteration        x        y

Initial Value     0         1

1          1         1

2          2         2

3          3         6

4          4        24

Since f (x, y) = xy and ∆ = 1, Euler’s method uses y(x + 1) = y + xy, where y(0) = 1.
The following shows the needed iterative steps.
• Iteration 1: We have x = 0 and y = 1 from initial value. Plugging this into the
above equation yields y(1) = 1 + 0 × 1 = 1. Thus, we have x = 1 and y = 1.
• Iteration 2: Since we have x = 1 and y = 1, plugging them into the above equation
yields y(2) = 1 + 1 × 2 = 2. Thus, we have x = 2 and y = 2.
• Iteration 3: Plugging x = 2 and y = 2 into the above equation yields y(3) =
2 + 2 × 4 = 6. Hence, we have x = 3 and y = 6.
• Iteration 4: Plugging x = 3 and y = 6 into the above equation yields y(4) =
6 + 3 × 6 = 24. Therefore, we have x = 4 and y = 24.

(b) [15 points] Solve the given equation on [0, 3] using Ralston’s second-order Runge-Kutta
method with ∆ = 1, and ﬁll your results in the following table.
Answer: The second order Runge-Kutta method uses the following formula:

y(x + ∆) = y(x) + ∆ [c1f (x, y) + c2f (x + p2∆, y + (a2,1f (x, y))∆)]

Ralston’s method is a special case in which c1 = 1 , c2 = 2 and p2 = a2,1 = 3 . Note that
3        3                 4
∆ = 1. The following table shows the computed results:

Iteration         x               y

Initial Value       0               1

3
1              1           2   = 1.5

81
2              2      16   = 5.0625

2025
3              3     64    = 31.640625
CS3911 Intro. to Numerical Methods with Fortran Final Solutions                                      7

If you wish to do hand calculation, Ralston’s formula can be simpliﬁed to reduce com-
putation eﬀort. Since f (x, y) = xy and ∆ = 1, Ralston’s formula can be rewritten as
follows:
1          2        3     3
y(x + 1) = y + f (x, y) + f x + , y + f (x, y)
3          3        4     4
The f (x, y) in 1 f (x, y) and 3 f (x, y) should be replaced by xy. This yields the following:
3              4

1     2    3     3
y(x + 1) = y + xy +   x + , y + xy
3     3    4     4
Since f (x, y) = xy (i.e., f (x, y) being the product of its two arguments), we have
3     3         3                 3
f x + , y + xy = x +                y + xy
4     4         4                 4
As a result, Ralston’s formula becomes
1     2    3               3
y(x + 1) = y + xy +   x+              y + xy
3     3    4               4
If you wish, the above can be simpliﬁed further, making calculations nearly eﬀortless.
Expanding the product yields:
2    3          3    25   1      1
x+         y + xy = xy + x2 y + y
3    4          4    24   2      2
Putting this back to the simpliﬁed Ralston’s formula yields
1       2      3    3
y(x + 1) = y + xy +        x + , y + xy
3       3      4    4
1       25     1      1
= y + xy + xy + x2 y + y
3       24     2      2
3      11       1 2
=   y + xy + x y
2       8       2
1         11
=   y 3 + x + x2
2          4
With this simpliﬁed form, calculations can be carried out easily.
• Iteration 1: Initially we have x = 0 and y = 1. Plugging x and y into the above
formula yields y(1) = 3 (i.e., x = 1 and y = 3 ).
2                     2
• Iteration 2: When iteration 2 starts, we have x = 1 and y = 3 . Plugging x = 1
2
and y = 3 into Ralston’s formula yields y(2) = 81 (i.e., x = 2 and y = 81 ).
2                                        16                   16
• Iteration 3: When iteration 3 starts, we have x = 2 and y = 81 . Plugging them
16
into Ralston’s formula gives y(3) = 2025 (i.e., x = 3 and y = 2025 ). This ﬁnishes
64                        64
the third iteration.

7. Random Numbers

(a) [10 points] Suppose we randomly draw a rectangle with side lengths a and b, both in
the range of [0, 1). What is the probability of having a rectangle whose area is larger
than of equal to 0.5?
The following shows a sequence of pre-generated random numbers. Each number has ﬁve
digits, and a decimal point is assumed to appear to the left of the ﬁrst digit. For example,
random number 53479 should be interpreted as 0.53479. Therefore, all random numbers
CS3911 Intro. to Numerical Methods with Fortran Final Solutions                                             8

are in the range of [0,1). Use Monte Carlo method with 12 runs to ﬁnd an approximation
of this probability. You should show all computation steps. Only providing an

x1 to x8      53479     81115    98036        12217        59526   40238     40577       39351
x9 to x16     43211     69255    97344        70328        58116   91964     26240       44643
x17 to x24    83287     97391    92823        77578        66023   38277     74523       71118

Answer: Since the lengths are random numbers in [0, 1), the given table of random
numbers can be used to estimate a probability.

Rectangle Number        Side a      Side b          Area = a × b     Area ≥ 0.5?
1             0.53479     0.81115            0.43379            No
2             0.98036     0.12217            0.11977            No
3             0.59526     0.40238            0.23952            No
4             0.40577     0.39351            0.15967            No
5             0.43211     0.69255            0.29926            No
6             0.97344     0.70328            0.68460           Yes
7             0.58116     0.91964            0.53446           Yes
8             0.26240     0.44643            0.11714            No
9             0.83287     0.97391            0.81114           Yes
10              0.92823     0.77578            0.72010           Yes
11              0.66023     0.38277            0.25272            No
12              0.74523     0.71118            0.52999           Yes

Of these 12 runs, ﬁve have an area larger than or equal to 0.5. Therefore, the estimated
probability is 5/12 = 0.41666 . . . ≈ 0.42. In fact, this estimated probability is a bit too
high. You may write a program to generate a large number of random rectangles and
the result should be approximately 0.15. If you know more about probability theory,
you should be able to determine this probability directly, because this is the probability
of two uniform random numbers. We will not provide any discussion in this direction as
this is not a topic in this course.

8. Eigenvalues and Eigenvectors

(a) [15 points] Use the power method to ﬁnd the largest eigenvalue and its corresponding
eigenvector of matrix A as shown below, and ﬁll the following table with your results.
The initial value (i.e., iteration 0) is z = [1, 1]T , and you only do two iterations (i.e.,
iterations 2 and 3).
           
1 2
           
A=               
4 3

Answer: This following table has the desired result:
CS3911 Intro. to Numerical Methods with Fortran Final Solutions                                                           9

Approx.                               Approx. Eigenvector

Iteration              Eigenvalue                                          x                 y

0                        1                                             1                 1

3
1                        7                            7    = 0.428571428                 1

33                                      17
2              7   = 4.714285714                      33   = 0.515151515                 1

Here is the details:
• Iteration 1: Since the initial value is z = [1, 1]T , we have w = A · z = [3, 7]T . The
maximum in w is 7, which is used to scale w. Thus, λ = 7 and z = w/7 = [ 3 , 1]T .
7
This completes iteration 1.
• Iteration 2: Since z = [ 3 , 1]T , we have w = A · z = [ 17 , 33 ]T . The maximum in w
7                                7 7
is 33 . We have λ = 33 and z = w/ 33 = [ 17 , 1]T . This completes iteration 2.
7               7                7    33

(b) [20 points] Use Jacobi method to ﬁnd all eigenvalues and their corresponding eigen-
vectors of the following symmetric matrix A:
                      √            √       
3             2            2
2            4       −    4
                                             
                                             
       √                                     
A=
           4
2           9
4
3
4


                                             
           √                                 
2       3            9
−       4        4            4
You should clearly provide all computation details, and match each eigen-
value with its corresponding eigenvector. Otherwise, you will risk low grade.
Additionally, you will receive zero point if you do not use Jacobi method.
One more note: Unless you are very good at hand calculation, use your cal-
culator to do this problem.
Answer: Since the maximum of oﬀ-diagonal entries is a2,3 = 3/4, we have
2a2,3      2× 3
tan(2θ) =                = 9 4 =∞
9
a3,3 − a2,2   4 − 4
and θ = π/4. The rotation matrix of iteration 1 is
                                                                                  
1          0                      0                        1       0        0
                                                                                     
                                                                                     
                                                                       √        √    
R1 =  0          cos(π/4)           sin(π/4)  =                   0          2        2   
                                                                       2        2    
                                                                                     
           √        √    
0 − sin(π/4) cos(π/4)                                                   2    2
0 −         2    2
The rotated A is computed as follows:
   3        1       
2        2   0
                                
                                
                                
A =    RT
1   · A · R1 =                  1
2
3
2   0 
                                
                                
0        0 3
CS3911 Intro. to Numerical Methods with Fortran Final Solutions                                                                10

Since A is not diagonal, it should be rotated again. In the second iteration, since the
maximum oﬀ-diagonal entry is a1,2 = 1 , we have
2

2a1,2      2× 1
tan(2θ) =               = 1 2 =∞
1
a2,2 − a1,1   2 − 2

Again, θ = π/4 and the second rotation matrix is
                                                     √       √                
2           2
cos(π/4)    sin(π/4) 0                              2           2        0
                                                                              
                                                                              
                                                     √         √              
R2 =  − sin(π/4) cos(π/4) 0  =                            − 2            2
0 
                                                       2           2          
                                                                              
                          
0            0               1                   0           0        1

The rotated A is computed as follows:
                       
1 0 0
                              
                              
                              
A = R T · A · R2 =  0
2                                       2 0 
                              
                              
0 0 3

Since A is diagonal, eigenvalues of matrix A are 1, 2 and 3.
The eigenvectors of A is computed as
                                                                                               √    √         
2    2
1      0              0                cos(π/4)                sin(π/4) 0                       2    2    0
                                                                                                               
                                                                                                               
                                                                                                          √    
R1·R2 =  0      cos(π/4)     sin(π/4) · − sin(π/4) cos(π/4) 0  =                                      −1     1     2   
                                                                                              2      2    2    
                                                                                                               
             √    
0 − sin(π/4) cos(π/4)                      0                       0           1                1
−1    2
2     2   2

The following table shows eigenvalues and their corresponding eigenvectors:

Eigenvalues            1                       2                   3
   √                     √                          
2                       2                      0
   2
              2
                   
                                           √        
 1                                                  
Eigenvectors        −2                          1              2       
                             2             2        
                                                    
   √        
1
2               −1
2                      2
2

More precisely, the corresponding eigenvector of the i-th eigenvalue is the i-th column
of R1 · R2 .

```
To top