Code is displaying graph but not plotting any points

조회 수: 2 (최근 30일)
Christopher Arreola
Christopher Arreola 2021년 12월 15일
댓글: Christopher Arreola 2021년 12월 15일
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)];

채택된 답변

KSSV
KSSV 2021년 12월 15일
편집: KSSV 2021년 12월 15일
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')

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Graph and Network Algorithms에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by