Problem 1: Let f : [a, b] → R, for a < b, be a bounded function and consider a point x and some spacing h > 0. We define the forward divided difference of f on x with spacing h by
δh,+f (x) = (f (x + h)h - f (x))/h (1)
provided that x + h ∈ [a,b].
As seen on Problem Sheet 3, the following results holds, if f is sufficiently smooth.
|δh,+ f(x) - f'(x)| ≤ h/2 maxξ∈ [a,b]|f''(ξ)|. (2)
The results of the above indicate that, as h is reduced, the finite difference approximation becomes more accurate. However, this results assumes exact arithmetic, we investigate what happens in a computer implementation.
(a) For f (x) = cos(x), write a piece of MATLAB code to evaluate the forward difference approximation δh,+ f(Π/4) for a given h. For h ∈ [10-12,10-2], plot the error |f' (Π/4) - δh,+ f (Π/4)| against h. What do you notice?
[Hint: Use logspace to uniformly divide the interval [10-12,10-2] with 1000 points on a logarithmic scale, then use loglog to plot the error against h, but plot the points individually, i.e., do not join the points with lines.]
(b) Show that, when using floating point arithmetic (assuming no underflow or overflow), the following holds:
(δh,+ f(x))^ := fl((fl(fl( f (x + h)) - fl( f (x)))/h) ≈ (f (x + h) - f (x))/h + (f (x + h)(∈1 + ∈3 + ∈4) - f(X)(∈2 + ∈3 + ∈4))/h
where |∈il < ∈M/2, for i = 1,.....,4 and ∈m <<1 is the machine epsilon.
[Hint: Use the results on Slide 30 of Section 3 and assume that x + h is stored exactly on the machine. Also, notice that ∈i∈j < ∈M, for i, j = 1, ... ,4.]
(c) Hence, show that, when floating point numbers are used, the truncation and round-off errors combine so that
|f'(x) - (δh,+ f(x))^| <≈ h/2 maxξ∈[x, x+h]|f''(ξ)| + 6(∈Mmaxξ∈[x, x+h]|f(ξ)|)/h
(d) For f (x) = cos(x), plot the theoretical error from part (c) against h on a pair of axes with logarithmic scales with x = Π/4. Compare your results with those from part (a).
[Hint: MATLAB's machine epsilon can be found by typing eps in the command window.]
Problem 2 We wish to solve the initial-value problem: For some interval [a, b] ⊂ R, find a function y : I -+ Rn such that
y'(x) = f (x, y(x)) and y(a) = y0 (3)
for some known vector function f : [a, b] x y([a, b]) → Rn and some values y0 ∈ Rn.
Consider a subdivision of the interval [a, b] into subintervals [xi-1, xi] of equal width, with
a = x0 < x1 < x2.... < xN = b
where xi = xi-1+ h, and h = (b - a)/N.
(a) Download newtonraphson.m, explicit_euler m and implicit_euler m from Blackboard. As written, the Euler functions will compute approimxate solutions to (3) with n = 1.
Apply the explicit and implicit Euler methods with f (x, y) = -1000y, y(0) = 1 and [a, b] = [0, 1].
(i) At what values of N (and hence h), does the explicit Euler method become unstable?
(ii) Is the implicit Euler method stable at these values? Does it ever become unstable?
(b) Apply the explicit and implicit Euler methods to the differential equation with f (x, y) = cos(4Π(x + y2)), on the interval [a, b] = [0, 1] and again let y(0) = 1.
(i) Using N = 5000, compute a good approximation to y(1) using either Euler method. Call this approxi-mation (1).
(ii) By choosing N = 10, 20, 40, 80, 160, 320, find the error |y(xN)-y~(1)| when both the implicit and explicit Euler method is used. Plot the error against N on a logarithmic scale. What is the order of convergence of both methods in terms of N and in terms of h?
(c) The Lotka-Volterra equations are a model for a predator-prey system in mathematical ecology. Let x(t) be the population of the prey at time t and y(t) be the population of the predator at time t. The equations read:
dx/dt = αx - βxy,
dy/dt = δxy - γy,
for some constants α, β, γ, δ.
Convert the implicit Euler code so that it can now solve systems of equations of the form (3), with n ≥ 1.
Let α = 0.1, β = 0.02, γ = 0.04 and δ = 0.004, x(0) = 20, y(0) = 2 and use your implicit Euler method to compute approximations to x(t) and y(t) for t ∈ [0, 1000] using a grid with time-step δt = 1. Plot x and y against t on a set of axis.
[Note, partial credit will be given if the explicit Euler method is used instead.]
Problem 3 - In this problem we use the finite difference method to solve the following convection-diffusion-reaction problem:
-u'' + cu' + du = f, x ∈ (a, b), (4)
u(a) = A, (5)
u(b) = B., (6)
where c > 0 and d > 0 are constants, f : [a, b] -> R is a known function and A and B are boundary conditions.
(a) Write a MATLAB function to approximate u(x) using the finite difference method, it should take as input the interval [a, b] and c, d, A and B and the number of subintervals N, such that
a = xo 1 = xo + h ..... < XN-2 + h < XN-1+ h = XN = b.
with h = (b - a)/N
(b) Let a = 0, b = 1, c = 1, d = 1, A = 0, B = 0 and f (x) = 1. Find the exact solution to the equation with these parameters. For N = 10, 20, 40, 80, 160, 320 run your code and compute the error:
EN = √(1/N∑i=0N(u(xi) - ui)2)
where ui is the finite difference approximation to u(xi).
Plot the solution for N = 3200 and give a plot of E against N on a set of axes with logarithmic scales. What is the order of the convergence of the method with respect to N?
Problem 4: The Euler methods can also be used to find solutions to problems of the form (4)-(6), this technique is called the shooting method. First, we convert the second order differential equation into a system of first order equations. By letting v = u', we have
du/dx = V,
dv/dx = cv + du - f (x).
We have u(a) = A, but we do not have information about v(a), however, we know that we need u(b) = B. We pick v(a) and solve to find u(b), then modify v(a) until u(b) = B. This is like changing the angle at which a gun is fired until the target is hit, hence the name.
For the same problem as in Problem 3(a), implement the shooting method with your Euler method from Problem 2(c) and find an approximation to u. For N = 20, 40, 80,160, 320 run your code and compute the error:
EN = √(1/N∑i=0N(u(xi) - ui)2)
Plot the solution for N = 320 and give a plot of E against N on a set of axes with logarithmic scales. What is the order of the convergence of the method with respect to N?
[Hint: Finding the correct v(a) is simple in this case because we have a linear differential equation, hence u(b) depends linearly on v(a).]