xuebaunion@vip.163.com

3551 Trousdale Rkwy, University Park, Los Angeles, CA

留学生论文指导和课程辅导

无忧GPA：https://www.essaygpa.com

工作时间：全年无休-早上8点到凌晨3点

微信客服：xiaoxionga100

微信客服：ITCS521

matlab代写-APMA0160

时间：2021-04-25

APMA0160 (A.Yew) Spring 2011

Numerical integration: Gaussian quadrature rules

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

arguments. The method underlying quadl is a “Gaussian quadrature rule”.

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.

学霸联盟

Numerical integration: Gaussian quadrature rules

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

arguments. The method underlying quadl is a “Gaussian quadrature rule”.

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.

学霸联盟