function retval = Wilkinson(N) % WILKINSON Compares lu-decomposition and qr-decomposition % This function compares lu-decomposition and qr-decomposition, % using the pathological matrix given by wilkinson % % usage: Wilkinson(n) and n is the maximal number of unknowns % for the system of equations to be solved % % Example: % % Wilkinson(50) % % Within the m-File, you can set the flags 'ITER' and 'DISPLAY' % for iterative refinement of the lu-decomposition and for % display of the residual norm, respectively % % (c) 1998, Rolf Krause, krause@math.fu-berlin.de if(nargin ~= 1), error('wrong number of arguments'); end if(N < 1), error('N < 2'); end DISPLAY = 0; % Display norm of residual Yes (1) / No (0) ITER = 0; % Iterative Refinment Yes (1) / No (0) lu_err = zeros(1,N); qr_err = zeros(1,N); cnd = zeros(1,N); cnd_r = zeros(1,N); for n=2:N n = fix(n); W = eye(n,n); W(:, n) = ones(n,1); for i = 1:n W(i, [1:i-1]) = - ones (1, i-1); end % assemble the rhs b = 1/n*linspace(0,n,n); b=b'; b_norm = norm(b, 2); [L,U,P] = lu(W); x = P*b; x = L\x; x = U\x; % iterative refinment if(ITER) r = b - W*x; r = L\r; r = U\r; x = x + r; end %x = x + W\(b-W*x); lu_err(n) = norm(W*x-b, 2) / b_norm; if(DISPLAY) disp(sprintf('dimension of linear system n = %d', n)); disp(sprintf('error of lu-decomposition: ||b-W*x||/||b|| = %12.7g', ... lu_err(n))); end [L, U, P] = qr(W); x = b; x = L'*x; x = U\x; x = P * x; qr_err(n) = norm(W*x-b, 2) / b_norm; if(DISPLAY) disp(sprintf('error of qr-decomposition: ||b-W*x||/||b|| = %12.7g\n', ... qr_err(n))); end [L, U, P] = lu(W); cnd(n) = cond(W); cnd_r(n) = cond(U); end currFig = gcf; clf; subplot(2,1,1); semilogy(lu_err([2:N]), 'r'); hold on; semilogy(qr_err([2:N]), 'g'); %semilogy(2*eps*[1:N].^3.*2.^[1:N], 'c'); legend('lu-decomposition', 'qr-decomposition',2);%, 'theoretical bound'); ax = axis; ax(1) = 2; axis(ax); title('relative error in the residual'); subplot(2,1,2); semilogy(eps*cnd(2:N)); title('eps*condition number of W'); figure(currFig+1); clf; semilogy(cnd(2:N), '+-'); hold on; semilogy(cnd_r(2:N), 'd-r'); legend('\kappa(W)', '\kappa(R)'); figure(currFig);