Code is displaying graph but not plotting any points

3 visualizaciones (últimos 30 días)
So this code is came from a textbook and was written in the late 90's, so syntax is not working properly in newer version of matlab. I need help figuring why none of the graphs are plotting any points.
function pivot
nn = [(5:5:100)];
iter = 5*ones(size(nn));
type = 1;
macheps=eps/2;
gpp = 0;
bndpp1 = 0;
bndpp2 = 0;
bndpp3 = 0;
bndpp4 = 0;
bndpp5 = 0;
bndpp6 = 0;
bndpp7 = 0;
bndpp8 = 0;
bndpp9 = 0;
bndpp10= 0;
bndpp11= 0;
errpp = 0;
normrpp= 0;
errppitref = 0;
normrppitref = 0;
gcp = 0;
bndcp1 = 0;
bndcp2 = 0;
bndcp3 = 0;
bndcp4 = 0;
bndcp5 = 0;
bndcp6 = 0;
bndcp7 = 0;
bndcp8 = 0;
bndcp9 = 0;
bndcp10= 0;
bndcp11= 0;
errcp = 0;
normrcp= 0;
errcpitref = 0;
normrcpitref = 0;
cnd1 = 0;
cnd2 = 0;
cnd3 = 0;
cnd4 = 0;
dim = 0;
j=0;
k=0;
for n = nn
k = k+1;
for i=1:iter(k)
j = j+1;
if (type ==1 )
a = randn(n,n);
end
if (type == 2 )
a = eye(n) + .0000001*randn(n,n);
dd = diag((10*ones(1,n)).^(14*(0:n-1)/(n-1)));
a = dd*a;
end
x = randn(n,1);
b = a*x;
[lpp,upp,ppp]=lu(a);
xpp = upp\(lpp\(ppp*b));
rppt = a*xpp-b;
dpp = upp\(lpp\(ppp*rppt));
xppitref = xpp - dpp;
gpp(j) = max(max(abs(upp)))/max(max(abs(a)));
bndpp1(j) = 3*n^3*gpp(j)*macheps;
bndpp2(j) = 3*n*norm(abs(lpp)*abs(upp),'inf')/norm(a,'inf')*macheps;
bndpp3(j) = norm(a*xpp-b,'inf')/(norm(a,'inf')*norm(xpp,'inf'));
bndpp5(j) = max(abs(a*xpp-b)./(abs(a)*abs(xpp)+abs(b)));
bndpp11(j)= max(abs(a*xppitref-b)./(abs(a)*abs(xppitref)+abs(b)));
errpp(j) = norm(xpp-x,'inf')/norm(xpp,'inf');
errppitref(j) = norm(xppitref-x,'inf')/norm(xpp,'inf');
[pcprow,lcp,ucp,pcpcol,err]=gecp(a);
xcp = pcpcol*(ucp\(lcp\(pcprow'*b)));
rcpitreft = a*xcp-b;
dcp = pcpcol*(ucp\(lcp\(pcprow'*rcpitreft)));
xcpitref = xcp - dcp;
gcp(j) = max(max(abs(ucp)))/max(max(abs(a)));
bndcp1(j) = 3*n^3*gcp(j)*macheps;
bndcp2(j) = 3*n*norm(abs(lcp)*abs(ucp),'inf')/norm(a,'inf')*macheps;
bndcp3(j) = norm(a*xcp-b,'inf')/(norm(a,'inf')*norm(xcp,'inf'));
bndcp5(j) = max(abs(a*xcp-b)./(abs(a)*abs(xcp)+abs(b)));
bndcp11(j)= max(abs(a*xcpitref-b)./(abs(a)*abs(xcpitref)+abs(b)));
errcp(j) = norm(xcp-x,'inf')/norm(xcp,'inf');
errcpitref(j) = norm(xcpitref-x,'inf')/norm(xcp,'inf');
inverse_a = pcpcol*(ucp\(lcp\(pcprow'*eye(n))));
cnd1(j) = norm(a,'inf')*norm(inverse_a,'inf');
[cnd2(j), v] = condest(a');
cnd3(j) = norm(abs(inverse_a)*(abs(a)*abs(xpp)+abs(b)),'inf') ...
/norm(xpp,'inf');
cnd4(j) = norm(abs(inverse_a)*(abs(a)*abs(xcp)+abs(b)),'inf') ...
/norm(xcp,'inf');
rpp = a*xpp - b;
normrpp(j)= norm(rpp,'inf');
bndpp4(j) = bndpp3(j) * cnd2(j);
bndpp6(j) = bndpp5(j) * cnd3(j);
bndpp7(j) = norm(abs(inverse_a)*abs(rpp),'inf')/norm(xpp,'inf');
bndpp8(j) = norm(abs(inverse_a)*(abs(rpp)+ ...
(n+1)*eps*(abs(a)*abs(xpp)+abs(b))),'inf')/ ...
norm(xpp,'inf');
rppitref = a*xppitref - b;
normrppitref(j)= norm(rppitref,'inf');
bndpp9(j) = norm(abs(inverse_a)*abs(rppitref),'inf')/norm(xpp,'inf');
bndpp10(j)= norm(abs(inverse_a)*(abs(rppitref)+ ...
(n+1)*eps*(abs(a)*abs(xppitref)+abs(b))),'inf')/ ...
norm(xpp,'inf');
rcp = a*xcp-b;
normrcp(j)= norm(rcp,'inf');
bndcp4(j) = bndcp3(j) * cnd2(j);
bndcp6(j) = bndcp5(j) * cnd3(j);
bndcp7(j) = norm(abs(inverse_a)*abs(rcp),'inf')/norm(xcp,'inf');
bndcp8(j) = norm(abs(inverse_a)*(abs(rcp)+ ...
(n+1)*eps*(abs(a)*abs(xcp)+abs(b))),'inf')/ ...
norm(xcp,'inf');
rcpitref = a*xcpitref - b;
normrcpitref(j)= norm(rcpitref,'inf');
bndcp9(j) = norm(abs(inverse_a)*abs(rcpitref),'inf')/norm(xcp,'inf');
bndcp10(j)= norm(abs(inverse_a)*(abs(rcpitref)+ ...
(n+1)*eps*(abs(a)*abs(xcpitref)+abs(b))),'inf')/ ...
norm(xcp,'inf');
dim(j) = n;
end
disp(['just finished n = ',int2str(n)])
end
disp('Figure 2.1 in textbook, when type=1')
figure(1)
plot(dim,gpp,'ow',dim,gcp,'+w');
hold on
grid
title('Pivot Growth Factors, GEPP = o, GECP = +')
xlabel('Matrix Dimension')
hold off
disp('Figure 2.2 in textbook, when type = 1')
hold on
figure(2)
subplot(211)
semilogy(dim,max(bndpp1,1e-18),'xw')
axis([0,max(nn),1e-18,1e-8])
semilogy(dim,max(bndpp2,1e-18),'+w')
semilogy(dim,max(bndpp3,1e-18),'ow')
grid
plot([0,100],macheps*[1,1],'-w')
title('Backward error in Gaussian elimination with partial pivoting')
xlabel('Matrix dimension')
text(min(nn),1e-9,'macheps = 2^(-53) = 1.1e-16')
subplot(212)
semilogy(dim,max(bndcp1,1e-18),'xw')
axis([0,max(nn),1e-18,1e-8])
semilogy(dim,max(bndcp2,1e-18),'+w')
semilogy(dim,max(bndcp3,1e-18),'ow')
plot([0,100],macheps*[1,1],'-w')
grid
title('Backward error in Gaussian elimination with complete pivoting')
xlabel('Matrix dimension')
text(min(nn),1e-9,'macheps = 2^(-53) = 1.1e-16')
hold off
figure(3)
semilogy(dim,cnd1,'+w')
grid
title('Condition number of A')
xlabel('Matrix dimension')
figure(4)
plot(dim,cnd2./cnd1,'+w')
hold on
grid
title('Ratio of Hager''s Estimate to True Condition Number')
xlabel('Matrix dimension')
disp('Figure 2.3 (type=1) or 2.4(a) (type=2), in textbook')
lmaxis = floor(log10( ...
max([bndpp4,bndcp4,bndpp6,bndcp6,bndpp7,bndcp7, ...
bndpp8,bndcp8,bndpp9,bndcp9,bndpp10,bndcp10])));
maxis = 10^lmaxis;
hold off
figure(5)
loglog(errpp,bndpp4,'ow')
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcp,bndcp4,'+w')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
for i=1:lmaxis+15
loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error vs. Error Bound (2.13), o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
disp('Figure 2.4(c) when type=2, in textbook')
hold off
figure(6)
semilogy(dim,bndpp5,'ow',dim,bndcp5,'+w'), grid
title('Componentwise relative backward error, o = GEPP, + = GECP')
xlabel('Matrix dimension')
hold off
figure(7)
loglog(errpp,bndpp6,'ow')
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcp,bndcp6,'+w')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
for i=1:lmaxis+15
loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error vs. Error Bound (2.7), o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
disp('Figure 2.4(b) when type=2 , in textbook')
hold off
figure(8)
loglog(errpp,bndpp7,'ow')
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcp,bndcp7,'+w')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
for i=1:lmaxis+15
loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error vs. Error Bound (2.14), o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
hold off
figure(9)
loglog(errpp,bndpp8,'ow'), grid
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcp,bndcp8,'+w'), grid
loglog([1e-16,maxis],[1e-16,maxis],'w-')
for i=1:lmaxis+15
loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error vs. Error Bound (2.14), |r| increased, o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
hold off
figure(10)
disp('Example 2.5 in textbook, when type=2')
loglog(errppitref,bndpp9,'ow')
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcpitref,bndcp9,'+w')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
for i=1:lmaxis+15
loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error after Iter Ref vs. Error Bound (2.14), o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
hold off
figure(11)
disp('Example 2.5 in textbook, when type=2')
loglog(errppitref,bndpp10,'ow')
axis([1e-16,maxis,1e-16,maxis]);
axis('square');
hold on;
loglog(errcpitref,bndcp10,'+w')
loglog([1e-16,maxis],[1e-16,maxis],'w-')
for i=1:lmaxis+15
loglog([1e-16,maxis/10^i],[1e-16*10^i,maxis],'--w')
end
title('True Error after Iter Ref vs. Error Bound (2.14), |r| increased, o = GEPP, + = GECP')
xlabel('True Error'), ylabel('Error Bound')
hold off
figure(12)
disp('Example 2.5 in textbook, when type=2')
semilogy(dim,bndpp11,'ow',dim,bndcp11,'+w'), grid
title('Componentwise backward error after Iter Ref, o = GEPP, + = GECP')
xlabel('Matrix dimension')
function [pr,l,u,pc,err]=gecp(a)
n=max(size(a));
pr=eye(n);pc=eye(n);
aa=a;
for i=1:n-1
am=max(max(abs(aa(i:n,i:n))));
[I,J]=find(abs(aa(i:n,i:n)) == am);
imax=I(1)+i-1;
jmax=J(1)+i-1;
if (imax ~= i)
temp = pr(:,i);
pr(:,i) = pr(:,imax);
pr(:,imax) = temp;
temp = aa(i,:);
aa(i,:) = aa(imax,:);
aa(imax,:) = temp;
end
if (jmax ~= i)
temp = pc(:,i);
pc(:,i) = pc(:,jmax);
pc(:,jmax) = temp;
temp = aa(:,i);
aa(:,i) = aa(:,jmax);
aa(:,jmax) = temp;
end
aa(i+1:n,i) = aa(i+1:n,i)/aa(i,i);
aa(i+1:n,i+1:n) = aa(i+1:n,i+1:n) - aa(i+1:n,i)*aa(i,i+1:n);
end
l=eye(n);l=l+tril(aa,-1);
u=triu(aa,0);
err=[norm(pr*l*u*pc'-a,1)/norm(a,1); cond(a); cond(l); cond(u)];

Respuesta aceptada

KSSV
KSSV el 15 de Dic. de 2021
Editada: KSSV el 15 de Dic. de 2021
plot([0,100],macheps*[1,1],'-w')
It seems you are using white color for plotting. White background doesn't show up the white color. Change it to other color.
Or set the background of figure to other color.
set(gcf,'color','k')

Más respuestas (0)

Categorías

Más información sobre Graph and Network Algorithms en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by