18.330 Problem set 5

Exercise 1: Finite differences via interpolation

Consider the simplest forward finite-difference approximation for ′( ):

( + ) − ( )

( ) =

ℎ

When we calculate this numerically, there are two sources of error: truncation error, coming from approximating the exact Taylor expansion with a finite piece of it, and floating-point roundoff error.

1. Suppose that we perturb the input, ℎ, by Δ . Calculate (analytically) an approximation to the (absolute) error Δ on the output to first order in Δ ; you should find that it grows like ℎ−1.

2. Suppose that the input perturbation size is mach; the error from [1] is then the roundoff error. Find an estimate for the value of ℎ at which the truncation error balances with the roundoff error, and find the size of the error there. Compare this with the plot that we did in class.

3. Consider an interval [ , ] and let be the midpoint of the interval. Use Lagrange interpolation to find an analytical expression for the unique

quadratic function that passes through ( , ( )), ( , ( )) and

( , ( )).

4. Use your result from [3] to derive the centered difference approximation for the derivative ′( ) in terms of equally-spaced points separated by a distance ℎ.

5. What approximation does it give for the second derivative ″( )?

6. Use [3] to find a backward difference expression for ′( ) using informa-tion at nodes −2 and −1.

7. Find numerically the rate of convergence of the results from [3] and [4] for equally-spaced points separated by a distance ℎ for the function sin(2 ) at = /4, for values of ℎ between 10−6 and 10−1.

Exercise 2: Integration using Simpson’s rule

In this problem we will derive the second-order Newton–Cotes quadrature rule,

known as Simpson’s rule, for calculating ∫ ( ) .

1

Suppose you are given an -point quadrature rule with nodes ( ) =0 and weights ( ) =0 for integrating over the interval $[-1, 1]. That is, the are + 1 points with −1 ≤ ≤ 1, and the are given to you such that

◦

• ( ) ≃ ∑ ( )

−1 =0

1. Construct a new quadrature rule for integrating over a general interval [ , ]. I.e., find ′ and ′ such that

∫

( ) ≃ ∑ ′ ( ′ )

=0

∫−1

( )

,

as follows:

Derive the basic second-order Newton–Cotes quadrature rule for

1

2. Use your results from [Exercise 1] to find the degree-2 polynomial 2 that agrees with at the three points = −1, 0, 1. (Leave your result in terms of the values (−1), (0) and (1).)

3. Integrate 2 interval [−1, 1] to approximate ∫ in terms of (−1, (0) and (1). Express this result as a quadrature rule.

4. Combine your answers to [2] and [3] to write down the basic (not compos-ite) Simpson’s rule for integrating f over [ , ].

5. Given an interval [ , ], subdivide it into equal-width subintervals, ap-ply the basic Simpson’s rule to integrate over each subinterval, and sum the results to obtain the composite Simpson rule for integrating over [ , ]. How many samples of f does this rule require? (Be careful not to overcount).

Exercise 3: Using Newton–Cotes methods

1. Implement the composite 0th (rectangle), 1st (trapezoid), and 2nd-order (Simpson) Newton–Cotes quadrature rules for integrating an arbitrary function over an arbitrary interval with + 1 points. Each should be a single function like rectangle(f, a, b, N).

Note that in the case of Simpson’s rule, we are using a total of + 1 points; how many intervals does this correspond to?

1Basic quadrature rules are for nodes distributed over a single interval; composite quadrature rules are obtained by splitting up a large interval into subintervals and using a basic rule on each subinterval.

2

2. Calculate ∫1 exp(2 ) using each method. Plot the relative error −1

( ) = approx( ) − exact

exact

as a function of for in the range [10, 106] (or use a higher or lower upper bound depending on the computing power you have available).

Do these errors correspond with the expectations from the arguments in lec-tures?

3. Do the same for ∫−21 exp(− 2) . Use the erf function from the Spe-cialFunctions.jl package to calculate the “exact” result. [Hint: Check carefully the help for that function to make sure of the definition used.]

4. We showed that the trapezium rule has error at most ( 2). Consider the following integral of a smooth, periodic function:

2

= ∫ exp(cos( ))

0

Plot the error in the trapezium rule in this case. How fast does it decay with ? [This will be important later in the course.]

Note that this integral can be calculated exactly as 2 0 (1), where 0 is a modified Bessel function, which can be evaluated at 1 using the SpecialFunctions.jl package as besseli(0, 1).

Exercise 4: Euler method for ODEs

1. Implement the Euler method in a function euler(f, x0, δt, t_final), assuming that 0 = 0. Your code should work equally well if you put vectors in, to solve the equation ẋ= f(x).

2. Use your code to integrate the differential equation

from

to

=

5

with initial condition

0=0.5

.

Plot the exact solution and

the analytical solution.

=

0.01, 0.05, 0.1, 0.5

On a

the numerical solution for values of

̇ = 2

. = 0

different plot show the relative error as a function of time, compared to

3. Do the same for ̇ = −2 with initial condition 0 = 3.

4. For the above two cases, calculate the error at = 5 when the time interval is split into pieces for 10 1000. Plot the error as a function of . What is the rate of convergence as ℎ → 0?

3

̈

A pendulum satisfies the ODE + sin( ) = 0, where is the angle with the vertical.

( , )̇ = 2

̇−

( )

Show analytically that the quantity (“energy”)

cos

5.

[ ( ( ), ( )) = 0]

is conserved along a trajectory, i.e. that

1

2

, so that

̇

0

̇0

))

.

̇

( ( ), ( )) = ( (

), (

6. Solve this equation using the Euler method for initial conditions (0, 1) to show that the energy is not conserved.

7. Draw the phase plane. Explain graphically what is happening in terms of what each step does.

8. Plot as a function of time for different values of . How fast does it grow? Explain this in terms of what happens at each step.

4