octave:2> format compact octave:3> octave:3> octave:3> A = hilb(12); octave:4> octave:4> octave:4> octave:4> A(1:4,1:4) ans = 1.00000 0.50000 0.33333 0.25000 0.50000 0.33333 0.25000 0.20000 0.33333 0.25000 0.20000 0.16667 0.25000 0.20000 0.16667 0.14286 octave:5> octave:5> octave:5> octave:5> n =12; octave:6> x = ones(n,1); octave:7> b = A*x; octave:8> A\b warning: matrix singular to machine precision, rcond = 2.49977e-17 warning: matrix singular to machine precision, rcond = 2.63277e-17 ans = 1.00000 1.00000 1.00000 0.99996 1.00024 0.99927 1.00123 0.99909 0.99966 1.00119 0.99915 1.00021 octave:9> octave:9> octave:9> b = b+1.e-12*randn(n,1); octave:10> A\b ans = 1.00001 0.99893 1.02670 0.71857 2.52855 -3.63271 8.74137 -4.65517 -1.16439 8.39105 -4.25810 2.30520 octave:11> b = A*ones(n,1); octave:12> c = b+1.e-14*randn(n,1); octave:13> A\c ans = 1.00000 0.99995 1.00127 0.98678 1.07115 0.78589 1.35559 0.74211 0.89972 1.33753 0.76080 1.05922 octave:14> cond(A) ans = 1.6769e+16 octave:15> cond(A,1) warning: matrix singular to machine precision, rcond = 2.63277e-17 warning: called from cond at line 75 column 12 ans = 3.7983e+16 octave:16> octave:16> octave:16> octave:16> n = 6; octave:17> A = vander([ones(6,1) .^2 ]); octave:18> A = vander([1:n] .^2 ]); parse error: syntax error >>> A = vander([1:n] .^2 ]); ^ octave:18> A = vander([1:n] .^2 )' A = Columns 1 through 5: 1 1024 59049 1048576 9765625 1 256 6561 65536 390625 1 64 729 4096 15625 1 16 81 256 625 1 4 9 16 25 1 1 1 1 1 Column 6: 60466176 1679616 46656 1296 36 1 octave:19> A = vander([1:n] .^2 ); octave:20> tau = 1.e-10; octave:21> octave:21> octave:21> octave:21> E = randn(n,n) .* abs(A); octave:22> octave:22> e = randn(n,1) .* abs(b); error: product: nonconformant arguments (op1 is 6x1, op2 is 12x1) octave:22> x = ones(n,1); octave:23> b = A*x; octave:24> e = randn(n,1) .* abs(b); octave:25> octave:25> B = A + tau*E; octave:26> c = b + tau*e; octave:27> octave:27> octave:27> y = B\c; octave:28> octave:28> %% exact forward error octave:28> octave:28> norm(x-y)/norm(x) ans = 0.0000079597 octave:29> octave:29> num = cond(A)*tau*(norm(e)/norm(b) +norm(E)/norm(A)) num = 0.022994 octave:30> den = 1-tau*norm(inv(A))*norm(E); octave:31> %% Theorem 2 bound octave:31> num/den ans = 0.023201 octave:32> %% backward error octave:32> norm(b-A*y)/(norm(E)*norm(y) + norm(e)) ans = 0.000000000063018 octave:33> diary off