Hello, if you have any need, please feel free to consult us, this is my wechat: wx91due
Numerical Analysis 1 – Week 1
Subjects covered
- Class introduction and mechanics.
- How to compute (evaluate) functions:
- Horner's method. Evaluating polynomials and series expansions.
- Wallis's' rule. Evaluating continued fractions.
- Computer representations of numbers: integers and floating point.
- Brief overview of computer internals.
- Computational error, stability and conditioning.
Readings
• “Testing math functions in Microsoft Cloud Numerics”, by Brorson, Moskowitz, and Edelman.
Linked on Canvas.
Problems
Problem 1
- Write a program to compute the value of ln(x) for real inputs 0 < x < 100 using the above series. Start with a stopping tolerance of 1e-6 to know when to terminate your sum. Your program should check the input, and depending upon the input value select either equation (1) or (2) to compute ln(x). Make sure you use Horner's rule for summing series – don't naively compute z^p at each iteration of your loop.
- Next write a test program which exercises your program by passing in values in the domain 1e- 2 < x < 100. The easiest way to generate such values is to use the Matlab function logspace(). Test the values you compute against those returned by the Matlab built-in function log(x). Use a relative testing tolerance to compare your results against Matlab's. I used reltol = 1e- 5*max(1,x). You may need to adjust the stopping tolerance in your ln(x) implementation in order to pass your test.
- Next add some tests to your test program which compute the round trip, exp(ln(x)), and verifies the return is equal to x (within testing tolerance). Test values in the domain 1e-2 < x < 100.
- Now modify your test program to test larger values of x, say x= 1e4. Please explain what happens in a README file accompanying your homework submission.
Problem 2
In problem 1 you used some series expansions to compute ln(x) for positive real x in the domain 0 < x < 100. When you tried using your implementation for larger x you should have found that your implementation either failed tests, or you needed to dramatically increase the number of terms in your sum to get good results. From the perspective of writing correct, performant numerical code, neither of these options is good.
The problem for large x is that series (2) converges slowly. One way to fix this situation is to use reduction on your input x prior to computing the series expansion. The following algorithm is similar to that used for the ln(x) implementation in fdlibm, which you can see at http://www.netlib.org/fdlibm/e_log.c. It works by dividing the input by 2 while accumulating an additive factor of ln(2) for each division. Then, when the input x is small enough, a series sum for y = ln(x) is used. The return is then y + n*ln(2), where n is the number of divisions. Here is the algorithm:
1. Accept input x and sanity test it (i.e. make sure it's real and positive).2. Initialize p2 = 0, lg2 = ln(2).3. Loop:
1. if x>1.8 then set x = x/2, and set p2 = p2+lg2.2. else break out of loop and go to step 7.
4. End of loop.5. Use the series expansion valid for 0 < x < 2 to compute y = ln(x).6. Return y + p2.
Problem 3
cf. the DLMF, https://dlmf.nist.gov/4.9 , equation 4.9.1.
Please write a program to compute the log of an input z using this continued fraction expression. Use Wallis’s algorithm to implement the series expansion. Stop your iteration when the difference between two successive convergents is less than 1e-8. Feel free to use my program computing sqrt(x) as a starting point.
Test your implementation using both comparison against Matlab and using the round-trip method. Test using the relative tolerance reltol = 1e-5*max(1,z).. This time test input values in the domain 1e-2 < z < 100.
Problem 4
