matlab代写-APMA0160

APMA0160 (A.Yew) Spring 2011
Matlab’s built-in numerical integration function [Q,fcount]=quad(f,a,b,tol) is essentially our
simp_compextr code with some further efficiency-enhancing features.
Note that quad requires scalar functions to be defined with elementwise operations, so f(x) = 2
1+x2
should be entered as f=inline(’2./(1+x.^2)’, ’x’) or f=@(x) 2./(1+x.^2)
The default tolerance for quad is 10−6.
Matlab has another efficient integration command called quadl, with the same input and output
Recall that each Newton–Cotes quadrature rule came from integrating the Lagrange polynomial that
interpolates the integrand f at n equally spaced nodes in the interval [a, b]. Thus, in general, we expect
the degree of exactness of the rule to be n− 1 (though, as we’ve seen, some rules turn out to have a
higher-than-expected degree of exactness ).
Is it possible to find quadrature rules that use n nodes but have degree of exactness higher than n−1?
Any (non-composite) quadrature rule can be written in the form
(?)
n∑
i=1
wi f(xi)
where the xi are nodes and the wi are non-zero constants called “quadrature weights”.
If the n nodes are required to be equally spaced, then once you are told whether or not to include the
endpoints (i.e. whether the rule should be closed or open), the positions of the nodes are all fixed, so
it only remains to find the n weights wi. These are given by
wi =
∫ b
a
Li(x) dx where Li(x) =
n∏
k=1, k 6=i
x− xi
xk − xi
But what if the nodes are not required to be equally spaced? If we put no restrictions on the positions of
the nodes, other than that they must lie in the interval [a, b], then (?) would contain 2n free parameters:
w1, . . . , wn and x1, . . . , xn. In general, 2n parameters could potentially define a polynomial of degree
up to 2n−1. Therefore, the highest degree of exactness we can expect to achieve with any quadrature
rule is 2n− 1.
To illustrate the procedure of finding such a quadrature rule with degree of exactness 2n − 1, let us
consider how to choose the wi and xi when n = 2 and the interval of integration is [−1, 1]. In other
words, we want to determine w1, w2 and x1, x2 such that
w1f(x1) + w2f(x2) gives the exact value of
∫ 1
−1
f(x) dx
whenever f is a polynomial of degree 2(2)− 1 = 3 or less.
It can be proved using the theory of orthogonal polynomials that, on the integration interval [−1, 1],
to achieve the maximum degree of exactness 2n− 1 with n nodes and n weights we should choose the
nodes x1, . . . , xn to be the roots (zeros) of the degree-n Legendre polynomial Pn(x); the weights are
then given by wi =
∫ 1
−1
n∏
k=1, k 6=i
x− xi
xk − xi dx, and so
n∑
i=1
wi f(xi) is an approximation of
∫ 1
−1
f(x) dx.
The Legendre polynomials can be defined via the recursive relation
Pk+1(x) =
2k + 1
k + 1
xPk(x)− k
k + 1
Pk−1(x) with P0(x) = 1, P1(x) = x.
The first few are P2(x) = 32(x
2 − 13), P3(x) = 52(x3 − 35x) and P4(x) = 18(35x4 − 30x2 + 3).
Each Pk(x) is a polynomial of degree k, and has k roots that all lie in the interval (−1, 1).
So, to find the quadrature rule with maximum degree of exactness using n nodes and n weights, in
principle we need to:
• Find the Legendre polynomial Pn(x).
• Find the roots of Pn(x) in (−1, 1); these will be our nodes x1, . . . , xn (and so the quadrature
rule will be an “open” rule).
• Find the Lagrange polynomial that interpolates the integrand f(x) at x1, . . . , xn.
• Integrate the Lagrange polynomial to determine the weights w1, . . . , wn.
Fortunately, the roots of the Legendre polynomials and their corresponding weights have been exten-
sively tabulated, so we can simply use these tables without redoing the calculations. The only issue
to be careful about is that the tabulated values were computed by taking [−1, 1] as the integration
interval, whereas in any given problem the integration interval may not be [−1, 1]. Therefore we need
to transform the tabulated values to analogous values on a general interval [a, b]. This is done through
the following change of variable:
−1 ≤ x˜ ≤ 1 ←→ a ≤ x ≤ b
if
x˜− (−1)
1− (−1) =
x− a
b− a
so that x =
b− a
2
x˜+
a+ b
2
∫ b
a
f(x) dx =
∫ 1
−1
f
(
b− a
2
x˜+
a+ b
2
)
dx
dx˜
dx˜
=
∫ 1
−1
f
(
b− a
2
x˜+
a+ b
2
)
b− a
2
dx˜
=
b− a
2
∫ 1
−1
f
(
b− a
2
x˜+
a+ b
2
)
dx˜
≈ b− a
2
n∑
i=1
wi f
(
b− a
2
x˜i +
a+ b
2
)
If we define h = b − a (the length of the interval) and c = a+b2 (the midpoint of the interval), then
the roots x˜i in [−1, 1] are transformed to the nodes xi in [a, b] via xi = h2 x˜i + c , and the quadrature
formula for approximating
∫ b
a f(x) dx will be
h
2 times the formula for approximating the equivalent
integral over [−1, 1].
The quadrature rules defined above, using the roots of Legendre polynomials as their nodes, are called
Gauss–Legendre rules. They have degree of exactness 2n − 1 (and order 2n). Gauss–Legendre
rules are open rules, and because the nodes are often positioned at irrational points in the interval,
when we code the adaptive composite rules by repeatedly halving the interval, many extra function
evaluations may need to be performed.
In order to make use of previously computed function values more efficiently, one could try to define a
closed Gaussian quadrature rule. Such a rule would have x1 = a and xn = b, and it turns out that the
appropriate choice of the n− 2 interior nodes should be the (transformed) roots of P ′n−1(x) in (−1, 1).
These roots and their associated weights are also available in tables, and the same transformation as
above would convert them to the corresponding nodes and quadrature formula on [a, b]. The rules
defined in this way are called Gauss–Lobatto rules. Note that because x1 and xn are now fixed,
the formula
∑n
i=1wif(xi) has only 2n− 2 free parameters, so the degree of exactness attainable by a
Gauss–Lobatto rule is 2n− 3 (order 2n− 2).
Gauss–Legendre nodes and coefficients
n Roots of Pn(x) Coefficients ci
2 −1/√3 1
1/

3 1
3 −√3/5 5/9
0 8/9√
3/5 5/9
4 −
√(
15 + 2

30
)/
35
(
18−√30)/36

√(
15− 2√30)/35 (18 +√30)/36√(
15− 2√30)/35 (18 +√30)/36√(
15 + 2

30
)/
35
(
18−√30)/36
Gauss–Lobatto nodes and coefficients
n x1, roots of P ′n−1(x), xn Coefficients ci
2 −1 1
1 1
3 −1 1/3
0 4/3
1 1/3
4 −1 1/6
−1/√5 5/6
1/

5 5/6
1 1/6
5 −1 1/10
−√3/7 49/90
0 32/45√
3/7 49/90
1 1/10
Matlab’s quadl employs a composite 4-node Gauss–Lobatto rule with some further efficiency-enhancing
features.  