%Nystrom Method
%brad nelson 2/4/12
%math 126 hw 4

%gauss scheme=0; composite trap=1
scheme=0;

%functions for p.3
ker=@(t,s) exp(t.*s);
f=@(x) exp(x)+(exp(x+1)-1)./(1+x);
fact=@(x) exp(x);

%functions for p.2
% ker=@(t,s) t.*s.^3;
% f=@(x) x.^0;
% fact=@(t) 1-5*t/24;

if scheme==0
nmax=20; %will test schemes using up to nmax nodes
gap=1;
else
nmax=1000;
gap=100;
end

e=zeros(1,numel(1:gap:nmax));
counter=0;
for n=1:gap:nmax %we will use n+1 points
counter=counter+1;

if scheme==0
%gauss nodes, weights
[xj,w]=gauss_nodes(0,1,n);
else
%composite trapezoid nodes, weights
[xj,w]=ct_nodes(0,1,n);
end


[x y]=meshgrid(xj,xj); % matricies with values to compute k matrix


k=ker(y,x); %starts k matrix
for i=1:n+1
    k(:,i)=k(:,i).*w(i); %weights k matrix for quadrature
end
b=f(xj); %solution vector
id=eye(size(k)); %identity matrix

uval=(id+k)\b; %solves for u values at nodes
uact=fact(xj); %actual values at nodes

%checks error at nodes

fn=f(xj); %starts reconstructed function
for i=1:n+1
   fn=fn-ker(xj,xj(i))*w(i)*uval(i);
end


e(counter)=max(abs(fact(xj)-fn));

%following produces plot of error as fn of n for 5 gauss nodes
if n==4 && scheme==0
    xs=0:0.01:1; %x values to check error at
    fn=f(xs); %starts reconstructed function
        for i=1:n+1
            fn=fn-ker(xs,xj(i))*w(i)*uval(i);
        end
    figure;
    hold on;
    plot(xs,abs(fact(xs)-fn));
    plot(xj,zeros(1,numel(xj)),'x');
    vals=[0.7 0.84 0.7 1.3 7.5]*10^-10;
    for i=1:5
       plot(xs,zeros(1,numel(xs))+vals(i),'.');
    end
    title('error for gaussian quad, n=5');
    xlabel('x');
    ylabel('error');
    hold off;
    print -depsc 3c.eps;

end

end

if scheme==1
figure;
loglog(1:gap:nmax,e);
xlabel('n');
ylabel('max error');

print -depsc 3a_ct.eps;
end

if scheme==0
figure;
semilogy(1:gap:nmax,e);
xlabel('n');
ylabel('max error');

print -depsc 3a_gauss.eps;
end