function errorResidualConditioningStabilityDemo(A)
% given matrix A, choose x randomly, define b = A*x and solve Ax = b
% by GEPP, giving xhat = A\b, then compute
% the relative residual norm ||r||/||b|| where the residual r = b-A*xhat,
Copyright By PowCoder代写 加微信 powcoder
% the relative error norm ||e||/||x||, where the error e = x-xhat,
% and the condition number of A, ||A||*||A^{-1}||
% the important point is that a stable algorithm gives a small relative
% residual norm, meaning not much more than the machine precision,
% but the relative error norm may be much bigger, by a factor of the
% condition number
% instructive to try this for
% (a) a random matrix
% (b) a matrix whose last column is a small perturbation of the sum of
% the first two columns — such a matrix is ill-conditioned since
% its columns are nearly linearly dependent
% (c) the incredibly rare example where GEPP is unstable, with n>=30
% (d) the same, with a tiny perturbation added: 1.e-10*randn(size(A))
[n,m] = size(A);
error(‘A must be square’)
% define solution x randomly and then define b = A*x
x = randn(n,1);
% solve Ax = b by GEPP giving computed solution xhat
xhat = A\b;
% compute residual and error
r = b-A*xhat; % this we can always compute in practice
e = x – xhat; % this we don’t usually know, but we’ve cooked it up here so we know it
% compute relative residual norm and relative error norm
rel_resid_norm = norm(r)/norm(b);
rel_error_norm = norm(e)/norm(x);
% compute condition number
condA = norm(A)*norm(inv(A)); % we usually don’t want to or cannot compute this
fprintf(‘rel resid norm is %g, rel error norm is %g, cond(A) is %g\n’, …
rel_resid_norm, rel_error_norm, condA);
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com