This is a sample solution to the Lab 18 problem of converting a Poisson solver into a Helmholtz solver.
From http://faculty.washington.edu/rjl/classes/am583s2014/notes/labs/lab18.html, R. J. LeVeque
In [ ]:
%pylab inline
In [ ]:
from scipy import sparse # to define sparse matrices
from scipy.sparse.linalg import spsolve # to solve sparse systems
The boundary value problem $u''(x) + k^2 u(x) = -f(x)$ is discretized to get the linear system of equations $U_{i-1} - 2U_i + U_{i+1} + k^2\Delta x^2 U_i = -\Delta x^2 f(x_i)$ for $i = 1, ~2,~\ldots,~n$. Recall that $U_0$ and $U_{n+1}$ are specified by the two boundary conditions required for this problem.
This function solves the BVP directly by setting up a tridiagonal matrix and using spsolve to efficiently solve the system.
In [ ]:
def solve_BVP_direct(f,x,u_left,u_right,k):
"""
On input:
f should be the function defining the ODE u''(x) = f(x),
x should be an array of equally spaced points where the
solution is to be computed,
u_left should be the boundary condition at the left edge x[0],
u_right should be the boundary condition at the right edge x[-1],
Returns:
u, vector of approximate values to solution at the points x.
"""
n = len(x) - 2 # number of interior points
print "Solving tridiagonal system with n =%4i" % n
dx = x[1]-x[0]
# Form the matrix:
d1 = ones(n)
d0 = (k**2 * dx**2 - 2) * ones(n)
A = sparse.spdiags([d1,d0,d1], [-1,0,1],n,n,format='csc')
# Form the right-hand side
fx = f(x)
rhs = -dx**2 * fx[1:-1]
rhs[0] = rhs[0] - u_left
rhs[-1] = rhs[-1] - u_right
# Solve the system for the interior points:
u_int = spsolve(A,rhs) # vector of length n-2
# Append the known boundary values:
u = hstack((u_left, u_int, u_right))
return u
Define a function to compare the computed solution to the true solution:
In [ ]:
def check_solution(x,u,u_true):
# u should be a vector of approximate solution values at the points x
# u_true should be a function defining the true solution.
plot(x,u,'ro')
plot(x,u_true(x),'b')
error = u - u_true(x)
error_max = abs(abs(error)).max()
print "Maximum error is %10.3e" % error_max
To test this out, define the function $f(x)$ and the true solution, and set the boundary values.
In [ ]:
k = 15.
u_left = 2.
u_right = -1.
c1 = u_left
c2 = (u_right - c1*cos(k)) / sin(k)
print "c1 = %g, c2 = %g" % (c1,c2)
f = lambda x: 0. * x
u_true = lambda x: c1 * cos(k*x) + c2 * sin(k*x)
n = 50 # number of interior points
x = linspace(0, 1, n+2)
u = solve_BVP_direct(f,x,u_left,u_right,k)
check_solution(x,u,u_true)
In [ ]:
nvals = 10*array([10,20,40,80,160])
errors = zeros(nvals.shape)
for j,n in enumerate(nvals):
x = linspace(0, 1, n+2)
u = solve_BVP_direct(f,x,u_left,u_right,k)
errors[j] = abs(u - u_true(x)).max() # maximum abs error over interval
print "\n n error"
for j,n in enumerate(nvals):
print "%4i %15.9f" % (n,errors[j])
loglog(nvals,errors,'o-')
Suppose we have 2 threads available and wish to parallelize the solution to the BVP. Solving a tridiagonal system by Gauss Elimination is inherently quite sequential, but we can try the following approach. Suppose we knew the value of the solution at the midpoint of the interval. Then we could solve 2 disjoint BVPs, one on the left half and one on the right half. In principle the two could be solved in parallel.
The problem is that we don't know the solution at the midpoint in general. But for this simple test problem we can try evaluating the exact solution at the midpoint to try this out...
In [ ]:
def solve_BVP_split(f,x,u_left,u_right,k):
n2 = len(x)
nmid = int(floor(n2/2.))
xhalf1 = x[:nmid+1]
xhalf2 = x[nmid:]
u_mid = u_true(x[nmid]) # Assumes we know true solution!!
print "u_mid = ",u_mid
uhalf1 = solve_BVP_direct(f,xhalf1,u_left,u_mid,k)
uhalf2 = solve_BVP_direct(f,xhalf2,u_mid,u_right,k)
u = hstack((uhalf1, uhalf2[1:]))
return u
In [ ]:
x = linspace(0, 1, 51)
u = solve_BVP_split(f,x,u_left,u_right,k)
check_solution(x,u,u_true)
Unfortunately, for this to work we need to know the solution in the middle. We do for this simple test problem but for a practical problem we would not. If you change the code defining solve_BVP_split above to set u_mid to some other arbitrary value, you will find that the solutions on each half do not match up smoothly and so the equation $u''(x) = -f(x)$ will not be satisfied at the midpoint.
Note that the finite difference equation at the midpoint is not being imposed. The equation $U_{i-1} - 2U_i + U_{i+1} = -\Delta x^2 f(x_i)$ is solved for $i=1,~2,~\ldots,nmid-1$ and for $i=nmid+1,~\ldots,~n$, but not at i = nmid.
This suggests an approach in which we guess the midpoint value, solve the two BVPs, and then determine the mismatch in the missing equation that we need to have satisfied at the midpoint. If this is non-zero then we can try to define a way to improve it, perhaps an iterative method, to drive the mismatch to zero.
First let's see how the mismatch varies as we change the guess for u_mid:
In [ ]:
def solve_BVP_split_mismatch(f,x,u_left,u_right,k):
n2 = len(x)
dx = x[1]-x[0]
nmid = int(floor(n2/2.))
xhalf1 = x[:nmid+1]
xhalf2 = x[nmid:]
x_mid = x[nmid]
u_mid_guesses = linspace(0,2,5)
mismatch = zeros(u_mid_guesses.shape)
for j,u_mid in enumerate(u_mid_guesses):
uhalf1 = solve_BVP_direct(f,xhalf1,u_left,u_mid,k)
uhalf2 = solve_BVP_direct(f,xhalf2,u_mid,u_right,k)
u = hstack((uhalf1, uhalf2[1:]))
plot(x,u)
mismatch[j] = (uhalf1[-2] + (k**2*dx**2 - 2.)*u_mid + uhalf2[1]) + dx**2 * f(x_mid)
print "With u_mid = %g, the mismatch is %g" % (u_mid, mismatch[j])
figure()
plot(u_mid_guesses, mismatch,'o-')
x = linspace(0, 1, 40)
solve_BVP_split_mismatch(f,x,u_left,u_right,k)
The top plot above shows the computed solution for various guesses at u_mid. The bottom plot shows the mismatch as a function of the guess u_mid.
Note that the behavior is linear. (You might try to use linear algebra to convince yourself that it should be!)
This means we don't need an iterative method to find the right value of u_mid. If we compute the solution for any two distinct values of u_mid and call the solutions obtained v0 and v1 and the mismatches G0 and G1, then we can solve for a value z such that z*G0 + (1-z)*G1 = 0 and the solution we desire is z*v0 + (1-z)*v1 by linearity. Note also that if the original boundary conditions are satisfied by both v0 and v1 then the convex combination z*v0 + (1-z)*v1 also satisfies these boundary conditions (and the ODE).
This suggests the following code, which solves each problem twice and the produces the right linear combination.
In [ ]:
def solve_BVP_split(f,x,u_left,u_right,k):
n2 = len(x)
dx = x[1]-x[0]
nmid = int(floor(n2/2.))
xhalf1 = x[:nmid+1]
xhalf2 = x[nmid:]
x_mid = x[nmid]
# solve the sub-problems twice with different values of u_mid
# Note that any two distinct values can be used.
u_mid = 0.
uhalf1 = solve_BVP_direct(f,xhalf1,u_left,u_mid,k)
uhalf2 = solve_BVP_direct(f,xhalf2,u_mid,u_right,k)
G0 = (uhalf1[-2] + (k**2*dx**2 - 2.)*u_mid + uhalf2[1]) + dx**2 * f(x_mid)
v0 = hstack((uhalf1, uhalf2[1:]))
u_mid = 1.
uhalf1 = solve_BVP_direct(f,xhalf1,u_left,u_mid,k)
uhalf2 = solve_BVP_direct(f,xhalf2,u_mid,u_right,k)
G1 = (uhalf1[-2] + (k**2*dx**2 - 2.)*u_mid + uhalf2[1]) + dx**2 * f(x_mid)
v1 = hstack((uhalf1, uhalf2[1:]))
z = G1 / (G1 - G0)
u = z*v0 + (1-z)*v1
return u
In [ ]:
x = linspace(0, 1, 80)
u = solve_BVP_split(f,x,u_left,u_right,k)
print x.shape, u.shape
check_solution(x,u,u_true)
The code above requires 4 calls to solve_BVP_direct and the four are independent and could be done in parallel. Suppose you have more threads or processes available. How would you use them? You could split the original interval into more pieces and have additional points where you would need to choose appropriate boundary conditions for each piece by requiring that the resulting solutions also satisfy the difference equations at these boundary points. It turns out this would give a linear system of equations to solve for these interface boundary values.
Another approach is to use recursion to split the interval into more pieces. Rather than calling solve_BVP_direct four times (which sets up and solves a tridiagonal system), we could recursively use the same strategy for each piece. This suggests the following recursive version:
In [ ]:
def solve_BVP_split_recursive(f,x,u_left,u_right,k):
n2 = len(x)
print "Entering solve_BVP_split_recursive with n2 =%4i, for x from %5.3f to %5.3f" \
% (n2, x[0], x[-1])
if n2 < 20:
# Stop recursing if the problem is small enough:
u = solve_BVP_direct(f,x,u_left,u_right,k)
return u
else:
# recursive
nmid = int(floor(n2/2.))
x_mid = x[nmid]
xhalf1 = x[:nmid+1]
xhalf2 = x[nmid:]
dx = x[1]-x[0]
u_mid = 0.
uhalf1 = solve_BVP_split_recursive(f,xhalf1,u_left,u_mid,k)
uhalf2 = solve_BVP_split_recursive(f,xhalf2,u_mid,u_right,k)
G0 = (uhalf1[-2] + (k**2*dx**2 - 2.)*u_mid + uhalf2[1]) + dx**2 * f(x_mid)
v0 = hstack((uhalf1, uhalf2[1:]))
u_mid = 1.
uhalf1 = solve_BVP_split_recursive(f,xhalf1,u_left,u_mid,k)
uhalf2 = solve_BVP_split_recursive(f,xhalf2,u_mid,u_right,k)
G1 = (uhalf1[-2] + (k**2*dx**2 - 2.)*u_mid + uhalf2[1]) + dx**2 * f(x_mid)
v1 = hstack((uhalf1, uhalf2[1:]))
z = G1 / (G1 - G0)
u = z*v0 + (1-z)*v1
return u
In [ ]:
x = linspace(0, 1, 40)
u = solve_BVP_split_recursive(f,x,u_left,u_right,k)
check_solution(x,u,u_true)
In [ ]:
In [ ]: