Question

Write a MATLAB code for the conjugate gradient method and apply it to solve the system Hx = b, where H is the n×n Hilbert matrix, and b is A times the vector of all ones, for (a) n = 4; (b) n = 8. Compare your numerical solutions with the exact solution (which is the vector of all ones), and report your numerical errors.

Answer 1

%%Matlab code for conjugate gradient method
clear all
close all

%Running the program with given parameters
%All A, x0, b and epsilon values
%Hilbert matrix of size (4X4)
A=hillbrt(4);
xx=ones(4,1);
b=A*xx;
x0=zeros(4,1);

fprintf('Example code for Conjugate gradient \n')
fprintf('The A matrix is \n')
disp(A)

fprintf('The b matrix is \n')
disp(b)

%Running the program with conjugate gradient function
[x,iter]=conjugategradient(A,x0,b);

fprintf('The solution matrix using conjugate gradient function is \n')
disp(x)

fprintf('Error in conjugate gradient method is %e.\n',norm(x-A\b))
  
%Running the program with given parameters
%All A, x0, b and epsilon values
%Hilbert matrix of size (4X4)
A=hillbrt(8);
xx=ones(8,1);
b=A*xx;
x0=zeros(8,1);

fprintf('Example code for Conjugate gradient \n')
fprintf('The A matrix is \n')
disp(A)

fprintf('The b matrix is \n')
disp(b)

%Running the program with conjugate gradient function
[x,iter]=conjugategradient(A,x0,b);

fprintf('The solution matrix using conjugate gradient function is \n')
disp(x)

fprintf('Error in conjugate gradient method is %e.\n',norm(x-A\b))  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Function for creating Hilbert Matrix
function A=hillbrt(n)
%creating Hilbert matrix of size nXn
    for i=1:n
        for j=1:n
            A(i,j)=1/(i+j-1);
        end
    end
end
  

%%Function for Conjugate Gradient method
function [x,iter]=conjugategradient(A,x0,b)
    epsi=10^-8;
    x=x0;
    r_old=b-A*x;
    p=r_old;
    epsi=epsi*norm(r_old);
    cnt=0;
    while norm(r_old)>epsi

        cnt=cnt+1;
        alpha=(r_old'*r_old)/(p'*A*p);
        x=x+alpha*p;
        r_new=r_old-alpha*A*p;
        beta=(r_new'*(r_new-r_old))/(r_old'*r_old);
        p=r_new+beta*p;
        r_old=r_new;

    end
    iter=cnt;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

       %%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%