%Q1

clear all
close all

rlist=0*[1:1:20];% empty lists to load r,p,q in
plist=rlist;
qlist=rlist;


rlist(1)=0.994;% n=0 values
plist(1)=1;
qlist(1)=1;

% we will write the LaTeX code for a table containing all enteries out to
% the file table.tex
fileID = fopen('table.tex','w');%open table.tex
fprintf(fileID,'%s \n','\begin{tabular}{ c | c|c|c|c}');%table header
fprintf(fileID,'%s \n','\hline');
fprintf(fileID,'%s \n','$n$ & $r_n$ & $p_n$ & $q_n$& $x_n$\\');
fprintf(fileID,'%s \n','\hline');
fprintf(fileID,'%s \n','0 & 0.994 & 1 & 1 & 1\\');% load n=0 values into table

plist(2)=0.497;%n=1 values for p&1
qlist(2)=0.497;

xlist=2.^(-[0:1:19]);% exact sequence for comparasion

for i=2:20% generate sequences
   rlist(i)=rlist(i-1)/2;% r_n
   
   
   if i>2 %already have p_1 & q_1
       plist(i)=3*plist(i-1)/2-plist(i-2)/2;% p_n
       qlist(i)=5*qlist(i-1)/2-qlist(i-2);% q_n       
   end
   
   % convert numbers to scientific notation decimal*10^(expo) for output
   expop=floor(log10(abs(plist(i))));% the power of 10 for p_i
   decimalp=plist(i)*10^(-expop);% the decimal number of 10 for p_i
   
   expoq=floor(log10(abs(qlist(i))));
   decimalq=qlist(i)*10^(-expoq);
   
   expor=floor(log10(abs(rlist(i))));
   decimalr=rlist(i)*10^(-expor);
   
   expox=floor(log10(abs(xlist(i))));
   decimalx=rlist(i)*10^(-expox);
   
   %output line of table
   lineoftex=strcat(num2str(i-1),...
       ' & $',num2str(decimalr),'\times 10^{',num2str(expor),'}$',...
       ' & $',num2str(decimalp),'\times 10^{',num2str(expop),'}$',...
       ' & $',num2str(decimalq),'\times 10^{',num2str(expoq),'}$',...
       ' & $',num2str(decimalx),'\times 10^{',num2str(expox),'}$\\');
   fprintf(fileID,'%s \n',lineoftex);

end

fprintf(fileID,'%s \n','\hline');%bottom of table
fprintf(fileID,'%s \n','\end{tabular}');

fclose(fileID);% close table.tex

xax=0:1:19;%xvalues for plot

figure(1)
%semilog plots of relative error
semilogy(xax,abs(xlist-rlist)./xlist,'o:b','LineWidth',2)
hold on
semilogy(xax,abs(xlist-plist)./xlist,'d-.r','LineWidth',2)
semilogy(xax,abs(xlist-qlist)./xlist,'^--m','LineWidth',2)
xlabel('$n$','interpreter','latex','FontSize', 18)
ylabel('$|$error$|/x_n$','interpreter','latex', 'FontSize', 18)

%inset plot of absolut error
axes('Position',[.25 .65 .3 .3])%inset plot possition
box on
semilogy(xax,abs(xlist-rlist),'o:b','LineWidth',2)
hold on
semilogy(xax,abs(xlist-plist),'d-.r','LineWidth',2)
semilogy(xax,abs(xlist-qlist),'^--m','LineWidth',2)
xlabel('$n$','interpreter','latex','FontSize', 14)
ylabel('$|$error$|$','interpreter','latex', 'FontSize', 14)

print( '-depsc2', 'Q1' )

%%
%%%%%%%%%%%%% Q2

%% changing the coeff of y
close all
A=[1/2,1/3;1/3,1/4];% coeff of (4)&(5)
B=[1/2,1/3;1/3,-1/2];% coeff of (6)&(7)
Am=[1/2,0.33;1/3,1/4];% modified coeff of (4)&(5)
Bm=[1/2,0.33;1/3,-1/2];% modified coeff of (4)&(5)
c=[1;-8];% right hand side of linear equations

x1exact=A\c%exact results
x2exact=B\c

x1witherror=Am\c%results with error
x2witherror=Bm\c

Rerror1=(x1exact-x1witherror)./x1exact% relative compont errors
Rerror2=(x2exact-x2witherror)./x2exact

R1=norm(x1exact-x1witherror)/norm(x1exact)%relative RMS error
R2=norm(x2exact-x2witherror)/norm(x2exact)

detA=det(A)%determinant of A
detAm=det(Am)

detB=det(B)%determinant of B
detBm=det(Bm)


figure(1)% plot of (4)&(5)

xA=[180:40:220];
plot(xA,3*(1-xA/2),'b-.','LineWidth',2)
hold on
plot(xA,(1-xA/2)/0.33,'r:','LineWidth',2)
plot(xA,-4*(8+xA/3),'k','LineWidth',2)
plot(x1exact(1),x1exact(2),'ok','LineWidth',2,'MarkerFaceColor','b')
plot(x1witherror(1),x1witherror(2),'dk','LineWidth',2,'MarkerFaceColor','r')
xlabel('$x$','interpreter','latex','FontSize', 18)
ylabel('$y$','interpreter','latex', 'FontSize', 18)
axis([185, 215,-320,-280])
set(gca,'XTick',[185 215])
set(gca,'YTick',[-320 -280])
print( '-depsc2', 'Q245a' )

figure(2)% plot of (4)&(5)

xB=[-6.1:0.3:-5.8];
plot(xB,3*(1-xB/2),'b-.','LineWidth',2)
hold on
plot(xB,(1-xB/2)/0.33,'r:','LineWidth',2)
plot(xB,2*(8+xB/3),'k','LineWidth',2)
plot(x2exact(1),x2exact(2),'ok','LineWidth',2,'MarkerFaceColor','b')
plot(x2witherror(1),x2witherror(2),'dk','LineWidth',2,'MarkerFaceColor','r')
xlabel('$x$','interpreter','latex','FontSize', 18)
ylabel('$y$','interpreter','latex', 'FontSize', 18)
axis([-6.1, -5.8,11.8,12.3])
set(gca,'XTick',[-6.1, -5.8])
set(gca,'YTick',[11.8,12.3])
print( '-depsc2', 'Q267a' )

%% changing the coeff of y
close all
A=[1/2,1/3;1/3,1/4];
B=[1/2,1/3;1/3,-1/2];
Am=[1/2,1/3;0.33,1/4];
Bm=[1/2,1/3;0.33,-1/2];
c=[1;-8];

x1exact=A\c
x2exact=B\c

x1witherror=Am\c
x2witherror=Bm\c

Rerror1=(x1exact-x1witherror)./x1exact
Rerror2=(x2exact-x2witherror)./x2exact

R1=norm(x1exact-x1witherror)/norm(x1exact)
R2=norm(x2exact-x2witherror)/norm(x2exact)

detA=det(A)
detAm=det(Am)

detB=det(B)
detBm=det(Bm)

