In [ ]: %matplotlib inline %precision 8
exam_1 March 2, 2018
from __future__ import print_function import numpy
import matplotlib.pyplot as plt
Before you turn this problem in, make sure everything runs as expected. First, restart the kernel (in the menubar, select Kernel → Restart) and then run all cells (in the menubar, select Cell → Run All).
Make sure you fill in any place that says YOUR CODE HERE or “YOUR ANSWER HERE”, as well as your name and collaborators below:
1 HW 3: BVP Problems II 1.1 Question 1
Consider the following BVP: where
and f (x) is defined by
uxx=f(x), x∈[−1,1] ux(−1) = 1 and u(1) = 0
0 x ∈ [−1, −a]
x/a+1 x∈[−a,0] f ( x ) = − x / a + 1 x ∈ [0, a]
0 x ∈ [a, 1]
(a) [10] Solve the BVP analytically.
(b) [20] Derive a method that solves the above BVP and write a function that accepts a as an input and computes the solution. You may use solve here but make sure to use at least 2nd order convergent approximations. Show that your method converges at the expected rate and handles variable values of a.
In [ ]: # INSERT CODE HERE
raise NotImplementedError(“Replace this statement with your solution.”)
(c) [5] Comment on the behavior with respect to varying a you see. 1
1.2 Question 2 – Higher Order Approximations
Let us consider a few approaches to achieving higher order convergence. For all of these questions consider the basic ODE
u′′(x)= f(x) x∈[0,L] u(0)=α, u(L)=β
For all of these you can solve the resulting system using an appropriate numpy or scipy command. (a) 4th order finite difference: [15] Using the following fourth order accurate finite difference
approximations to write a function that solves u′′ (x) = f (x).
u′′(xj) ≈ −Uj−1 + 16Uj−1 − 30Uj + 16Uj+1 − Uj+2
12∆x2
u′′(x1) ≈ 10U0 −15U1 −4U2 +14U3 −6U4 +U5
12∆x2
recalling that the forward and backward differences come in pairs with minor changes.
In [ ]: def high_order_differencing(m, f, alpha, beta, L): “””
Solve u''(x) = f(x) using 4th order differencing methods
:Input:
- *m* (int) Number of points to use in the approximation (not including boundaries
- *f* (func) The right hand side function in the ODE
- *alpha* (float) Left hand boundary condition u(0)
- *beta* (float) Right hand boundary condition u(L)
- *L* (float) Length of domain such that x \in [0, L]
:Output:
- *x* (ndarray) Points at which the solution is found AND the boundaries
- *U* (ndarray) Solution at the x points and the boundaries.
“””
# INSERT CODE HERE
raise NotImplementedError(“Replace this statement with your solution.”)
return x, U
In [ ]: alpha = -5.0
beta = 3.0
L = 3.0
u_true = lambda x, alpha, beta, L: numpy.exp(x) + (beta – alpha + 1.0 – numpy.exp(L)) / x, U = high_order_differencing(10, lambda x: numpy.exp(x), alpha, beta, L)
error = numpy.linalg.norm(U – u_true(x, alpha, beta, L), ord=numpy.infty)
print(error)
assert(error < 1e-3)
print(“Success!”)
2
(b) Deferred Corrections: [15] One way to obtain higher order results is to compute using a second order method but twice. The first step is to compute the solution to the system
AU = F.
Recall now that we know that the global error is defined as E = U − Uˆ and satisfies
We have an approximation to τ as
AE = −τ.
∆x2 (4) 4 τj = 12 u (xj)+O(∆x )
which depends on the exact solution but we have already computed an approximation to that and can use an approximation to the fourth derivative of U to compute τ. We can then solve the new system of equations and update the previous problem with the estimate of the global error we found.
In general deferred correction methods solve the original problem at lower order and then use this solution to compute an estimate of the truncation error τ, find the approximation of the global error E, and use this to update the solution.
For this question implement this method for u′′(x) = f(x) again. Note that there may be a significant short cut assuming that the 4th derivative of u(x) is smooth and using the original equation u′′(x) = f(x). Another hint, this shortcut is very similar to the one we found for the 9-point Laplacian. Do NOT assume that you know f (x) before hand.
In [ ]: def deferred_correction(m, f, alpha, beta, L): “””
Solve u''(x) = f(x) using a deferred correction method
:Input:
- *m* (int) Number of points to use in the approximation (not including boundaries
- *f* (func) The right hand side function in the ODE
- *alpha* (float) Left hand boundary condition u(0)
- *beta* (float) Right hand boundary condition u(L)
- *L* (float) Length of domain such that x \in [0, L]
:Output:
- *x* (ndarray) Points at which the solution is found AND the boundaries
- *U* (ndarray) Solution at the x points and the boundaries.
“””
# INSERT CODE HERE
raise NotImplementedError(“Replace this statement with your solution.”) return x_bc, U
In [ ]: alpha = -5.0
beta = 3.0
L = 3.0
3
u_true = lambda x, alpha, beta, L: numpy.exp(x) + (beta – alpha + 1.0 – numpy.exp(L)) / x, U = deferred_correction(10, lambda x: numpy.exp(x), alpha, beta, L)
error = numpy.linalg.norm(U – u_true(x, alpha, beta, L), ord=numpy.infty)
print(error)
assert(error < 1e-3) print(“Success!”)
(c) [5] Confirm the order of convergence for the two methods above. In [ ]: # INSERT CODE HERE
raise NotImplementedError(“Replace this statement with your solution.”) 1.3 Question 3
[30] Your colleague knocks on your door and presents you with a linear system of equations de- rived from a finite difference approximation of a PDE and asks whether they should expect Ja- cobi, Gauss-Seidel or SOR to converge for the resulting linear system. Describe how you would go about determining whether these methods would converge and at what rate you would expect them to do so. Use both analytical and computation arguments to justify your answer along with some test cases.
In [ ]: # INSERT CODE HERE
raise NotImplementedError(“Replace this statement with your solution.”)
4