Math 104A Computation Homework # 2
Jea-Hyun Park
Due: Nov 29, 11:59 PM
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.
1. composite trapezoidal rules (“Trapezoidal” below) (textbook p. 481)
2. composite Simpson’s rule (“Simpson” below) (textbook p. 484)
3. 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 C1– function, and (a scaled) normal distribution density.
f(x) = cos(πx), g(x) =
x
2
(x ≥ 0)
−x
2
(x < 0)
, h(x) = exp(
−x
2
2
). (1)
• (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.
1
– 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
r =
log(e1/e2)
log(N2/N1)
.
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:
#Subintervals Trapezoidal Simpson Gauss3
consine C1 function normal dist consine C1 function normal dist consine C1 function normal dist
1
2
4
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 Simpson Simpson Simpson Gauss3 Gauss3 Gauss3
Integrand consine C1 function normal dist consine C1 function normal dist consine C1 function normal dist
1
2
4
– 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:
2
#Subintervals consine C1 function normal dist
Trapezoidal Simpson Gauss3 Trapezoidal Simpson Gauss3 Trapezoidal Simpson Gauss3
1
2
4
Integrand consine consine consine C1 function C1 function C1 function normal dist normal dist normal dist
Method Trapezoidal Simpson Gauss3 Trapezoidal Simpson Gauss3 Trapezoidal Simpson Gauss3
1
2
4
• (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]
• This tells Wolfram Alpha that find R 1.1
−1
exp(−x
2
/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.
3
• (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.
4
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.
5
