AERO 300 Laboratory 6
Numerical Differentiation and Newton-Cotes Integration (Quadrature)
The pre-lab assignment (Section 4) is due at the beginning of this lab (can be hand-written or
typed). The lab (Section 5) will be due before your next lab section (.zip file submitted to
PolyLearn).
1 Objectives
This lab covers concepts and methods used for approximating the derivative of a function using
finite difference equations. In addition, this lab covers numerical techniques for integrating a
function with Newton-Cotes equations. In both cases, this lab assumes the closed form of the
function is unknown; the function is only known at discrete points. In this lab, you will:
• Develop code to approximate the first and second derivative of data representing the
displacement of a falling object.
• Use the Rectangle, Trapezoid and Simpson’s rule methods of quadrature to integrate
data from an accelerometer to calculate speed and displacement.
• In both cases, the data sets include measurement noise.
2 Introduction
Quadrature is the computation of a univariate definite integral. It can refer to either numerical
or analytic techniques; one must gather from context which is meant. The term refers to the
geometric origin of integration in determining the area of a plane figure by approximating it
with squares.
There are several reasons for ca
ying out numerical integration. First, the integrand f(x) may be
known only at certain points, such as obtained by sampling. Some embedded systems and
other computer applications may need numerical integration for this reason. Even if a formula
for the integrand may be known, it may be difficult or impossible to find an antiderivative which
is an elementary function. An example of such an integrand is ?(?) = ?−?
2
, the antiderivative
of which cannot be written in elementary form. For some functions, it may be possible to find
an antiderivative symbolically, but it may be easier to compute a numerical approximation than
to compute the antiderivative. That may be the case if the antiderivative is given as an infinite
series or product, or if its evaluation requires a special function which is not available.
Numerical Differentiation is the “inverse” process of numerical integration. All the reasons to
take a numerical derivative ca
y over from those for numerical integration. In addition, many
equations which govern the motion of fluids and structure are most easily solved numerically
using finite difference equations in multiple dimensions (space and time). We will return to this
use of finite difference equations later in the course. In the meantime, we will use finite
difference schemes to perform straight ahead estimate of the derivatives of data.
2.1 Numerical Differentiation
As discussed in class, numerical differentiation is accomplished by using the Taylor series expansion to
develop a finite difference equation. With respect to e
ors in the approximation of the derivative,
different finite difference formulas were also developed. Table 1 below summarizes those formulas.
Two-Point Forward-Difference
?′(?) ≈
?(? + ℎ) − ?(?)
ℎ
Three-Point Centered Difference
?′(?) ≈
?(? + ℎ) − ?(? − ℎ)
2ℎ
Three-Point Centered-Difference (2nd
derivative) ?
′′(?) ≈
?(? − ℎ) − 2?(?) + ?(? + ℎ)
ℎ2
Oftentimes you are given data and you want to take the derivative of the data. In this case, you may
know what the function should look like, but the data contains noise. In this lab, you will investigate
how the presence of noise impact the numerical differentiation of data.
2.2 Numerical Integration – Quadrature
As with numerical differentiation, in some instances, you have measurement data and you want to
calculate the integral of the data. Again, you may know what the function you are integrating looks like,
ut the data contains noise.
In addition to calculating the total integral using composite Newton-Cotes methods, we can compute
the integral at each time step, store that information so we can then plot the integral of a function
versus time. For example, using an accelerometer, you can measure the acceleration of a falling object
every 0.1 seconds. If you integrate the measurement, you would then have the speed of the object. If
you integrate again, you would then have the displacement. We call this “propagation” of the object. In
this case, we know from physics what the speed and displacement of the object should be as a function
of time.
We can use Newton-Cotes quadrature methods to propagate the speed and displacement by applying
any of the Newton-Cotes formulas to integrate the speed or displacement over one time step and add
the calculated integral to the previous value. For example, if the speed and displacement of an object
are 0
?
?
and 3 ? at ? = 0, and we measure the acceleration as ?(0) = 9.85
?
?2
and ?(0.1) = 9.75
?
?2
, we
can then estimate the speed and displacement using the trapezoid rule as:
?(0.1) = ?(0) +
ℎ
2
(?(0) + ?(0.1)) = 0 +
0.1
2
XXXXXXXXXX) = 0.98
?
?
?(0.1) = ?(0) +
ℎ
2
(?(0) + ?(0.1)) = 3 +
0.1
2
XXXXXXXXXX) = 0.98
?
?
= 3.049 ?
In general, we could continue this process for each time step we measure the acceleration and use the
following propagation equations:
??+1 = ?? +
ℎ
2
(?? + ??+1)
??+1 = ?? +
ℎ
2
(?? + ??+1)
In these equations, the ? represents the time step, and ℎ is the time between samples. As described
elow, we can use the same process using any of the quadrature methods.
Rectangle Rule
In integral calculus, the rectangle method uses an approximation to a definite integral, made by
finding the area of a series of rectangles. Either the left or right corners, or top middle of the
oxes lie on the graph of a function, with the bases run along the x-axis. Each of these schemes
is shown below in figure 1.
Figure 1. Left, midpoint, and right rectangular Integration schemes
To propagate a measurement using the left-rectangular rule, we can use the following formula:
?
?+1
= ?
?
+ ℎ?(?
?
)
(1)
where we have taken
∫ ?(?)??
??+1
??
≈ ℎ?(??)
Trapezoid Rule
The trapezoid rule uses trapezoids, as opposed to rectangles to approximate an integral. An
advantage of the trapezoid rule is that the sign of the e
or of the approximation is easily
known. An integral approximated with this rule on a concave-up function will be an
overestimate because the trapezoids include all the area under the curve and extend over it.
Using this method on a concave-down function yields an underestimate because area is
unaccounted for under the curve, but none is counted above. If the interval of the integral
eing approximated includes an inflection point, then the e
or is harder to identify.
Figure 2 – Trapezoid rule scheme
To propagate a measurement using the trapezoid rule, we can use the following formula:
??+1 = ?? +
ℎ
2
(?(??) + ?(??+1))
Simpson’s Rule
Simpson's rule is another method for approximating definite integrals.
It works in a similar way to the trapezoidal rule except that the
integrand is approximated to be a quadratic rather than a straight line
within each subinterval.
To propagate a measurement using Simpson’s rule, we can use the
following formula:
??+1 = ??−1 +
ℎ
3
(?(??+1) + 4?(??) + ?(??−1))
2.3 Accuracy
Each method has an accuracy associated with it. The rectangle rule, for example, is accurate for
a zeroth order polynomial. Using the approximation on any higher order polynomial results in
an e
or; if the rectangle rule was used to approximate a first order polynomial the e
or would
e on the order of h. In other words, the rectangle rule is a first order accurate method. A
similar argument can be made for the trapezoid and Simpson's rules.
3 Help
Potentially useful terms for this lab:
• for
• linspace
• ones
• plot
• grid
• hold
• legend
• title
• xlabel
• ylabel
4 Pre-lab Assignment
NOTE: This must be done before lab. If it is not completed, you will not be allowed into lab.
Figure 5: Simpson's approximation
Andy G Cristales
Use the concept of left-rectangular propagation from equation (1) to show that the total integral from ?
to ? is given as:
∫ ?(?)??
??=?
?1=?
≈ ℎ?(?) + ℎ ∑ ?(??)
?−1
?=2
5 Lab Assignment
Complete the following problems with a single script and custom functions as appropriate.
Important lines of code should include descriptive comments. You are encouraged to work in
groups, however the code you submit must be your own. All graphs/figures/sketches should
include a grid, labels, legend (if necessary), and the appropriate fontsize/linewidth/markersize.
1. Numerical Differentiation
a. Load the data file dispData.mat. The data comes from measuring the
displacement of a falling object versus time.
. Use the two-point forward-difference and the three-point centered-different
schemes to estimate the derivative of the given data.
c. Use the three-point centered-difference scheme to determine the 2nd derivative
of the data.
d. Assuming the initial displacement and speed of the object are zero meters and
zero meters per second, plot the following (plot the displacements, speeds and
accelerations on separate subplots):
i. The expected displacement, speed and acceleration versus time.
ii. The measured displacement and estimated speed and acceleration versus
time.
e. What is the effect of taking a numerical derivative of noisy data? Can you think of
a better method for estimating the acceleration of gravity from this data?
2. Numerical Integration – Quadrature
a. Load the data file accelData.mat. The data comes from measure tin the
acceleration of a falling object versus time.
. Use the rectangular, trapezoid and Simpson’s rules to integrate the data twice,
hence propagating the speed and displacement of the object.
c. Assuming the initial displacement and speed of the object are zero meters and
zero meters per second, plot the following (plot the