Added matrix conditioning demos

drreynolds · Jan 25, 2021 · 6d426ba · 6d426ba
1 parent 82753b4
commit 6d426ba
Show file tree

Hide file tree

Showing 2 changed files with 72 additions and 0 deletions.
diff --git a/IllConditioning/Demo.m b/IllConditioning/Demo.m
@@ -0,0 +1,29 @@
+% Script to demonstrate effects of matrix conditioning in floating-point arithmetic.
+%
+% Daniel R. Reynolds
+% SMU Mathematics
+% Math 4315
+clear
+
+% set matrix sizes for tests
+nvals = [6 8 10 12 14];
+
+% run tests for each matrix size
+for n = nvals
+
+  % create matrix, solution and right-hand side vector
+  A = hilb(n);
+  x = rand(n,1);
+  b = A*x;
+
+  % ouptut condition number
+  fprintf('Hilbert matrix of dimension %i: condition number = %g\n', n, cond(A));
+
+  % solve the linear system
+  S = warning('off','MATLAB:nearlySingularMatrix');
+  x_comp = A\b;
+
+  % output relative solution error
+  fprintf('  relative solution error = %g\n', norm(x-x_comp)/norm(x))
+
+end
diff --git a/IllConditioning/Demo.py b/IllConditioning/Demo.py
@@ -0,0 +1,43 @@
+#!/usr/bin/env python3
+#
+# Script to demonstrate effects of matrix conditioning in floating-point arithmetic.
+#
+# Daniel R. Reynolds
+# SMU Mathematics
+# Math 4315
+
+# imports
+import numpy
+
+# utility function to create Hilbert matrix of dimension n
+def Hilbert(n):
+  H = numpy.zeros([n,n])
+  for i in range(n):
+    for j in range(n):
+      H[i,j] = 1/(1+i+j)
+  return H
+
+
+# set matrix sizes for tests
+nvals = [6, 8, 10, 12, 14]
+
+# run tests for each matrix size
+for n in nvals:
+
+  # create matrix, solution and right-hand side vector
+  A = Hilbert(n)
+  x = numpy.random.rand(n,1)
+  b = A@x
+
+  # ouptut condition number
+  print("Hilbert matrix of dimension", n, ": condition number = ",
+        format(numpy.linalg.cond(A),'e'))
+
+  # solve the linear system
+  x_comp = numpy.linalg.solve(A,b)
+
+  # output relative solution error
+  print("  relative solution error = ",
+        numpy.linalg.norm(x-x_comp)/numpy.linalg.norm(x))
+
+# end of script