Labs for October 7th, 2015
Exercise 0.1. Implement all the methods presented in the lecture, i.e.:
numpy import dot, array, zeros_like, allclose, max, abs
numpy.linalg import norm
copy import deepcopy
scipy.linalg import lu, inv
1
from
from
from
from
6
def noPivot(A, row):
return row + 1
def partialPivot(A, row):
return abs(A[row:,row]).argmax() + row
11
def scaledPivot(A, row):
scales = max(abs(A), 1)
return abs(A[row:,row]/scales[row:]).argmax() + row
16
21
26
31
def GaussJordan(A, b, pivoting = noPivot):
(rows, cols) = A.shape
for row in range(0, rows-1):
pivot = pivoting(A, row)
if abs(A[pivot, row]) < 1e-8: raise ValueError()
if pivot != row:
A[[row, pivot],:] = A[[pivot, row],:]
b[[row, pivot]] = b[[pivot, row]]
for i in range(row+1, rows):
if abs(A[row, row]) < 1e-8: raise ValueError()
factor = A[i, row] / A[row, row]
A[i, row+1:rows] = A[i, row+1:rows] - factor*A[row, row+1:rows]
b[i] = b[i] - factor*b[row]
for k in range(rows-1,-1,-1):
b[k] = (b[k] - dot(A[k, k+1:rows], b[k+1:rows])) / A[k, k]
return b
def GaussJordanPartial(A, b):
return GaussJordan(A, b, pivoting = partialPivot)
36
def GaussJordanScaled(A, b):
return GaussJordan(A, b, pivoting = scaledPivot)
def LU(A, b):
L, U = lu(A, permute_l=True)
1
41
46
51
56
61
66
71
76
81
# Solve Ly = B for y ("forward")
y = zeros_like(b)
for m, b in enumerate(b.flatten()):
y[m] = b
if m:
for n in xrange(m):
y[m] -= y[n] * L[m, n]
y[m] /= L[m, m]
# Solve y = Ux for x ("backward")
x = zeros_like(b)
for midx in xrange(B.size):
m = b.size - 1 - midx # backwards index
x[m] = y[m]
if midx:
for nidx in xrange(midx):
n = b.size - 1 - nidx
x[m] -= x[n] * U[m, n]
x[m] /= U[m, m]
return x
def Jacobi(A, b, tol = 1e-10, limit = 100):
x = zeros_like(b)
for iteration in range(limit):
next = zeros_like(x)
for i in range(A.shape[0]):
s1 = dot(A[i, :i], x[:i])
s2 = dot(A[i, i + 1:], x[i + 1:])
next[i] = (b[i] - s1 - s2) / A[i, i]
if allclose(x, next, atol=tol): break
x = next
return x
def GaussSeidel(A, b, tol = 1e-10, limit = 100):
x = zeros_like(b)
for iteration in range(limit):
next = zeros_like(x)
for i in range(A.shape[0]):
s1 = dot(A[i, :i], next[:i])
s2 = dot(A[i, i + 1:], x[i + 1:])
next[i] = (b[i] - s1 - s2) / A[i, i]
if allclose(x, next, rtol=tol): break
x = next
return x
.
2
Test the methods on some well-known systems, trace the computation step by step, and try to explain
the results, as much as you can:
2
7
# An example from https://en.wikipedia.org/wiki/Gauss-Seidel_method
# with solution [-38 29]
A1 = array([[2.0, 3.0], [5.0, 7.0]])
b1 = array([11.0, 13.0])
A2 = array([[0.0, 1.0], [1.0, 1.0]])
b2 = array([1.0, 1.0])
# An example of Trefethen and Bau, with solution [-1, 1]
A3 = array([[1e-20, 1.0], [-1.0, 0.0]])
b3 = array([1.0, 1.0])
12
17
22
27
# Another example of Trefethen and Bau, with loss of n-bits of precision
# in n x n system. Try the same structure for larger n.
A4 = array([[ 1.0, 0.0, 0.0, 0.0, 1.0],
[-1.0, 1.0, 0.0, 0.0, 1.0],
[-1.0, -1.0, 1.0, 0.0, 1.0],
[-1.0, -1.0, -1.0, 1.0, 1.0],
[-1.0, -1.0, -1.0, -1.0, 1.0]])
b4 = array([0.0, 0.0, 0.0, 0.0, 0.0])
# An example from Faddeev and Faddeeva
A5 = array([
[5.0,
7.0,
6.0,
5.0],
[7.0,
10.0,
8.0,
7.0],
[6.0,
8.0,
10.0,
9.0],
[5.0,
7.0,
9.0,
10.0]])
b5 = array([23.0, 32.0, 33.0, 31.0])
.
Exercise 0.2. Consider the following routine for generating systems with an arbitrary condition number
c and prescribed solution (all ones):
2
7
import random
from numpy.linalg import svd
def randomIllConditioned(n = 10, c = 10000):
A = random.rand(n,n)
U, s, V = svd(A)
ss = zeros((n,n))
for j in range(n):
ss[j, j] = s[0] - j * (s[0] - s[0] / c) / (n - 1)
next = dot(dot(U, ss), V.T)
b = dot(next, ones(n))
return (next, b)
.
Plot the error of each method presented in class as a function of the condition number.
3
Exercise 0.3. Extend your experiments above with the conjugate gradient:
4
9
def CG(A, b, x0, tol = 1.0e-10, limit = 100):
x = x0
r0 = b - dot(A, x)
p = r0
for i in xrange(limit):
a = float(dot(r0.T, r0)/dot(dot(p.T, A), p))
if allclose(x, x + p*a, rtol=tol): break
x = x + p*a
r = r0 - dot(A*a, p)
b = float(dot(ri.T, ri)/dot(r0.T, r0))
p = ri + b * p
r0 = ri
return x
.
Additionally, plot the number of iterations and time needed to solve a random system for the three
iterative methods.
Hint: use module random, function random.rand(n,n).
4