Description

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 information 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 ∫ 𝑓(𝑥)𝑑𝑥.
𝑎
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. Construct a new quadrature rule for integrating over a general interval
[𝑢,𝑣]. I.e., find 𝑡′𝑘 and 𝑤𝑘′ such that
𝑣 𝑁
∫ 𝑓(𝑥)𝑑𝑥 ≃ ∑𝑤𝑘′𝑓(𝑡′𝑘)
𝑢 𝑘=0
Derive the basic second-order Newton–Cotes quadrature rule for 𝑓(𝑥)𝑑𝑥, as follows:
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 composite) Simpson’s rule for integrating f over [𝑢,𝑣].
5. Given an interval [𝑎,𝑏], subdivide it into 𝑁 equal-width subintervals, apply 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?
2. Calculate (2𝑥)𝑑𝑥 using each method. Plot the relative error
𝐼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 lectures?
3. Do the same for )𝑑𝑥. Use the erf function from the SpecialFunctions.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 𝑥̇ = 2𝑥 from 𝑡 = 0 to 𝑡 = 5 with initial condition 𝑥0 = 0.5. Plot the exact solution and the numerical solution for values of 𝛿𝑡 = 0.01,0.05,0.1,0.5. On a different plot show the relative error as a function of time, compared to the analytical solution.
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?
A pendulum satisfies the ODE 𝜃 +̈ sin(𝜃) = 0, where 𝜃 is the angle with the vertical.
5. Show analytically that the quantity (“energy”) 𝐸(𝜃, ̇ ̇
is conserved along a trajectory, i.e. that 𝑑𝑡𝑑 [𝐸(𝜃(𝑡),𝜃(𝑡)) = 0]̇ , 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.