Code is displaying graph but not plotting any points
2 views (last 30 days)
Show older comments
Christopher Arreola
on 15 Dec 2021
Commented: Christopher Arreola
on 15 Dec 2021
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)];
0 Comments
Accepted Answer
More Answers (0)
See Also
Categories
Find more on Graph and Network Algorithms in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!