driver_invheat.m

nonlinear eq routines in matlab

%
%  driver_invheat.m
% 
%  
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  March 29, 2001
%


clear all

% define A, the SVD U, D, V of A, the perturbed right hand side bdelta
% and the rank r of A as globals.

global A U D V bdelta


n     = 100;
kappa = 1;
[A,b,x] = heat(n,kappa);


% compute Singular Value Decomposition
[U,D,V] = svd(A);

figure(1); clf
semilogy(diag(D))
ylabel('Singular Values \sigma_i')

% set rank of A
r = size(A,2);

disp(' ')
disp(' hit return to continue ')
pause
clear delta


% add noise to the right hand side
delta  = 1.e-3;
bdelta = b.*( ones(n,1) + delta*rand(n,1));

disp([' solution of regularized LLS with noise levels delta= ',num2str(delta)])
disp('  ')


for i = 1:8
    % set regularization parameter lambda
    lambda(i) = 10^(-i);

    % compute regularized solution
    w = diag(D)./(diag(D).^2 + lambda(i)*ones(n,1));
    xlambda = V(:,1:r)*(w(1:r).*(U(:,1:r)'*bdelta));
    xerr(i) = norm(x-xlambda);
    res(i)  = norm(b-A*xlambda);  
end

figure(2)
loglog(lambda,xerr,'--'); hold on
loglog(lambda,res,'-');   hold off
xlabel('\lambda')
legend('|| x-x_{\lambda}||_2', '|| b^\delta-Ax_{\lambda}||_2')



disp(' Hit return to continue')
pause

f     = 'invheat';

fplot(f,[0,1]);
xlabel(' \lambda')
ylabel(' || b^\delta - A x_\lambda ||_2 ')


llow  = 0;
lup   = 1;

[llow, lup, lambda, ithist, iflag] = bisect( f, llow, lup );


disp(['The Bisection method returned with iflag = ',int2str(iflag)])
disp(' lambda estimated by the Bisection method:')
disp(['lambda= ', num2str(lambda)])