Math 104A Computation Homework # 2: Quadrature rules
Math 104A Computation Homework 2
Math 104A Computation Homework # 2
1 Main problem
We have studied some quadrature rules and their theoretical rates of convergence (or error terms). In this homework, you compute and report actual rates of convergence of some quadrature rules when they are implemented for selected integrands. Also, you interpret the result in your own language.
• (Quadrature rules) Your program computes the following three quadrature rules for integrals on [−1, 1.1] with respect to constant weight w(x) ≡ 1.
- composite trapezoidal rules (“Trapezoidal” below) (textbook p. 481)
- composite Simpson’s rule (“Simpson” below) (textbook p. 484)
- composite Gaussian quadrature with three nodes (“Gauss3” below), i.e., n = 2 so we have Gaussian nodes x0 , x1 , and x2 for each subinterval (textbook p. 496 provides a version on a single interval [−1, 1]. You will need to do some work to obtain a composite version from it. In particular, you will need to apply change of interval.)
• (Integrands) Implement these quadrature rules for the following three integrands: consine, a C 1 – function, and (a scaled) normal distribution density.
• (Experiment) Compute the errors caused by each quadrature rule for each integrand (i.e., the absolute value of the difference between approximation and the true value) implemented with a number of different resolutions. – Here, the resolution means the number of subintervals involved in composite rules. (2i for i = 0, 1, 2, · · · , 9 is a good option.) For example, if you decide to use 10 different resolutions, you will have total 3 × 3 × 10 errors from each combination of a method, an integrand, and a resolution. – To measure the error, you will need true values of the integrals. Think about how you can obtain those. If you have no idea, see the Section 2 for tips. • (Report) Report a table that contains the rate of convergence of each quadrature rule for each integrand as the resolution gets one step higher. – For example, if you decide to use 10 different resolutions, you will have 9 convergence rate for each combination of a quadrature rule and an integrand because you have 9 steps of refinement. If you are not sure what is going on, read ‘how to compute the rate of convergence’ below first, and possibly all other instructions, then come back here. – How to compute the rate of convergence? We will discuss this in class in a bit more detail. But the summary is the same. Suppose you have implemented a numerical method using two different resolutions N1 and N2 , where N2 is the higher resolution than N1 , say N2 = 2N1 . If you have obtained errors e1 and e2 using the resolutions N1 and N2 respectively, the rate of convergence based on these two results is given by See the class notes for details. • (Summary) data – Target integral information (interval [−1, 1.1]), Three integrands listed above parameters to choose – Resolutions (i.e., an array of the number of subintervals for composite rules; 2i for i = 0, 1, 2, · · · , 9 recommended) output – MethodIntegrand.csv: A csv file (comma separated values) that contains the rates of convergence so that one can easily compares them across integrands for each method used. Ideally, we want something like this:
However, csv files work the best when the data is complete and in the simplest rectangular form like the following. Thus, feel free to submit in this form: Method Trapezoidal Trapezoidal Trapezoidal – IntegrandMethod.csv: A csv file (comma separated values) that contains the rates of convergence so that one can easily compares them across methods for each integrand. For the same reason as above about csv files, feel free to report csv-friendly table: • (Submission): Please, submit your result as indicated below. Part of the credit is assigned to submitting things as instructed. – Source code: (.m, .py, or .ipynb; in particular, NO .pdf) – two output files: MethodIntegrand.csv, IntegrandMethod.csv – Interpretation: interpretation.txt (A text file that describes how you interpret the result) – (Optional): ExtraExperiment.txt – Interpretation on an additional experiment: Try the same experiment with for integrals from -1 to 1. (You must change the true integrals accordingly.) What do you observe? What is your interpretation?
2 Tips for this homework
2.1 Obtain true integrals
Although there should be many different ways, using Wolfram Alpha is one good way. Mathematica is case-sensitive. So be careful of capitalization. Also, use brackets and curly brackets correctly. • Type in the Wolfram Alpha: N[Integrate[Exp[-x*x/2], {x, -1, 1.1}], 20] R 1.1 • This tells Wolfram Alpha that find −1 exp(−x2 /2)dx up to 20 decimal digits • 20 decimal digits will do its job as the true value. • Do similar things for other integrands.
2.2 Tools for reporting
As I showed in class, it is very good to know data processing tools when you report the results. However, it takes time to learn those tools. Thus, consider learning and using them if you can budget some time. If you bring this up in office hours, I am happy to discuss it. • (Python) Use the pandas package (along with numpy). Among many great tools that belong to it, the most relevant things include: multiindex, pivot, stack, unstack, sort_values, to_csv, to_excel, etc. • (Matlab) Use table object. Note that the way Matlab approaches tables are quite different from how pandas treats data. I believe they are developed based on different philosophy. The most relevant tools include: splitvars, mergevars, stack, unstack, writetable, etc.
3 Tips for programming
3.1 Possible time savers
Here I share some issues/tricks you may want to know. Of course, your code may not experience any of these. But I am trying to save your time in case you decide to use the same ideas/tools and/or have the same issues.
Common • You can use for loop to conduct many different experiments in a single line. To do that, you can collect functions in a list (Python) or a cell array (Matlab). (See below for a bit more details for each language.) For example, if you use Python (Matlab), you can create a list (cell array) named fn of functions and a list (cell array) mthd of quadrature rules. Then, you can use for loop to conduct all experiments: “for i from 0 to 2 (from 1 to 3 for Matlab), and for j from 0 to 2 (from 1 to 3 for Matlab), compute the numerical integration of fn[j] using mthd[i].
Python • Depending on the version of Python, 5/3 may “equal” 1 instead of 1.66666 · · · . You can avoid this by typing “5./3”. • Depending on how you construct the resolutions, they may be considered float64 to Python rather than int (integer), which may cause error when Python carries out composite quadrature rules. Consider data type casting by NN = NN.astype(int), where NN is the name of the numpy array containing resolutions. • (List of functions) To collect functions, you simply gather them by creating a Python list of the function names. For example, you may want something like the following: fn = [f, g, h] mthd = [Trapez, Simpson, Gauss3]. If needed, Google “list of functions.”
Matlab • (cell array) You will want to use cell array for this homework. If you don’t know yet, study bare minimum of cell arrays. (Again, Google is a good friend). a quick piece of basic is that you have to use curly brackets to access its entry: e.g., fn{1}. if you use fn(1), Matlab will throw an error. • (cell array of functions) To collect functions, you have to collect so-called function handles of the Matlab functions. Anonymous functions’ names themselves are function handles (things you define using the syntax f = @(x) 2*x; ). But the functions you define using function output = f(x) .... end, you have to refer to it by @f. And you have to collect them in a cell array.
For example, if you have defined the integrands using anonymous functions and quadrature rules using regular functions, you may need something like the following: fn = {f, g, h} mthd = {@Trapez, @Simpson, @Gauss3}.
3.2 General comments
• Use Google to look up tools you need. For example, you may want to Google how to save numpy (Matlab) array as a csv file. This course does not give tutorials for everything, but encourages to independently find necessary programming tools. • General tips from Computational HW1 also apply. I don’t include them to avoid repetitions. For example, conducting unit tests and, if doable, vectorized coding is encouraged. (But if vectorizing makes things confusing, go for for loops first, then consider vectorizing it. Another example is is how to construct “anonymous functions.” Whenever you need, consider going over tips from Computational HW1.
