In [ ]: %matplotlib inline
%precision 8
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 two-dimensional Poisson problem defined as
∇2u = f(x,y) Ω = [0,1] [0,1]
with Dirichlet boundary conditions u(x,y)j¶Ω = 0.
Notethatthenotation¶ΩoftenreferstotheboundaryofΩandu(x,y)j¶Ω thesolutionevaluted
at the boundary.
(a) [5] If we wanted to consider solutions of the form
u(x,y) = (y3 y)(cos(2px) 1)
what should we require f(x,y) to be? Is this consistent with the boundary conditions?
(b)[10]Implementafinitedifferencemethodfortheproblemabovewithgridspacing∆x = ∆y
using a 9-point Laplacian.
In [ ]: # Suggested modules to use to construct matrix A
# You do not have to use these, they just may be helpful
import scipy.sparse as sparse
import scipy.sparse.linalg as linalg
def solve_9point(m, f):
# INSERT CODE HERE
raise NotImplementedError("Replace this statement with your solution.")
return x, y, U
1
In [ ]: f = lambda x, y: -(20.0 * y**3 + 9.0 * numpy.pi**2 * (y - y**5)) * numpy.sin(3.0 * numpy.pi * x)
u_true = lambda x, y: (y - y**5) * numpy.sin(3.0 * numpy.pi * x)
x, y, U = solve_9point(100, f)
X, Y = numpy.meshgrid(x, y)
error = numpy.linalg.norm((x[1] - x[0]) * (u_true(X, Y) - U), rd=1)
print(error)
assert error < 1e-3
print("Success!")
(c) [5] Show that the method is second-order accurate by preforming a convergence study (i.e.
plot the error vs. ∆x and compare this to the slopes for first and second order accurate methods).
In [ ]: # INSERT CODE HERE
raise NotImplementedError("Replace this statement with your solution.")
(d) [15] Show that the 9-point Laplacian can be written as a 5-point Laplacian plus a finite
difference approximation of 16h2uxxyy +O(h4).
(e) [10] Modify your function to use the trick introduced in class that will cause the 9-point
Laplacian stencil to become 4th order accurate. Show that this is true via a convergence study.
In [ ]: # INSERT CODE HERE
raise NotImplementedError("Replace this statement with your solution.")
1.2 Question 2
Let us consider some modifications to the Gauss-Seidel method we discussed and when one ver-
sion might be better than the other.
(a) [5] The Gauss-Seidel method for the discretization of u′′(x) = f(x) takes the form
if we assume we are marching forwards across the grid, for i = 1, 2, ..., m. We can also define a
backwards Gauss-Seidel method by setting
Show that this is a matrix splitting method of the type described in the lecture notes and what the
splitting is.
(b) [5]ImplementthebackwardsGauss-Seidelmethodandshowthattheconvergenceisatthe
samerateastheforwardversion(dothiscomputationally,youdonotneedtodothisanalytically).
Use the same expected iteration count as with the forwards method.
In [ ]: def solve_BGS(a, b, alpha, beta, m, f):
# INSERT CODE HERE
raise NotImplementedError("Replace this statement with your solution.")
return x_bc, U
In [ ]: # Problem setup
a = 0.0
b = 1.0
alpha = 0.0
beta = 3.0
f = lambda x: numpy.exp(x)
u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)
x, U = solve_BGS(a, b, alpha, beta, 150, f)
error = numpy.linalg.norm((x[1] - x[0]) * (u_true(x) - U), rd=1)
print(error)
assert(error < 1e-4)
print("Success!")
(c) [10] Modify the code so that it solves the steady-state problem
ϵu′′(x) u′(x) = f(x)
we had before with the same boundaries and f(x). This time however use the one-sided approxi-
mation to the first derivative U
Test both forward and backward Gauss-Seidel for the resulting linear system. Use the value
ϵ = 0.01, m = 100, and a maximum number of iterations of 250. Plot the convergence of both
approaches.
In [ ]: a = 0.0
b = 1.0
alpha = 1.0
beta = 3.0
epsilon = 0.01
u_true = lambda x: alpha + x + (beta - alpha - 1.0) * (numpy.exp(x / epsilon) - 1.0) / (numpy.exp(1.0 / epsilon) - 1.0)
# INSERT CODE HERE
raise NotImplementedError("Replace this statement with your solution.")
(d) [5] Explain intuitively why sweeping in one direction works so much better than in the
other.
Hint: Note that this equation is the steady equation for an advection-diffusion PDE
ut(x,t) + ux(x,t) = ϵuxx(x,t) f(x).
You might consider how the methods behave in the case ϵ = 0.
1.3 Question 3 - Poisson on a Rectangle
(a) [15] Write a function that solves the Poisson problem
∇2u = 54ex+y/2
3
on Ω = [0,1] [0,2] using a five point stencil. Use the true solution of the problem
u(x,y) = ex+y/2
to set the boundary conditions. Allow for ∆x and ∆y to be non-uniform.
In [ ]: import scipy.sparse as sparse
import scipy.sparse.linalg as linalg
def solve_poisson(X, Y, f):
# INSERT CODE HERE
raise NotImplementedError("Replace this statement with your solution.")
return U
In [ ]: # Grid including boundaries
# Note that the x and y arrays include the boundary location
m = 10
n = 25
delta_x = 1.0 / (m + 1)
delta_y = 2.0 / (n + 1)
x = numpy.linspace(0.0, 1.0, m + 2)
y = numpy.linspace(0.0, 2.0, n + 2)
X, Y = numpy.meshgrid(x, y)
# Transpose these so that the coordinates match up to (i, j)
X = X.transpose()
Y = Y.transpose()
f = lambda x, y: 1.25 * numpy.exp(x + y / 2.0)
u_true = lambda x, y: numpy.exp(x + y / 2.0)
U = solve_poisson(X, Y, f)
assert numpy.linalg.norm(max(delta_x, delta_y) * (u_true(X, Y) - U), rd=1) < 1e-3
print("Success!")
fig = plt.figure()
fig.set_figwidth(fig.get_figwidth())
axes = fig.add_subplot(1, 2, 1, aspect='equal')
plot = axes.pcolor(X, Y, U, vmax=7.0, vmin=0.0, cmap=plt.get_cmap("Blues"))
fig.colorbar(plot, label="$U$")
axes.set_title("Computed Solution")
axes.set_xlabel("x")
axes.set_ylabel("y")
axes = fig.add_subplot(1, 2, 2, aspect='equal')
plot = axes.pcolor(X, Y, u_true(X, Y), vmax=7.0, vmin=0.0, cmap=plt.get_cmap("Blues"))
fig.colorbar(plot, label="$u(x,t)$")
axes.set_title("True Solution")
axes.set_xlabel("x")
axes.set_ylabel("y")
plt.show()(b) [15] Keeping ∆x constant study the convergence behavior. as ∆y ! 0. Do the same for ∆x.
Comment on what you observe.
In [ ]: # INSERT CODE HERE
raise NotImplementedError("Replace this statement with your solution.")