Effacer les filtres
Effacer les filtres

Code is displaying graph but not plotting any points

2 vues (au cours des 30 derniers jours)
Christopher Arreola
Christopher Arreola le 15 Déc 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)];

Réponse acceptée

KSSV
KSSV le 15 Déc 2021
Modifié(e) : KSSV le 15 Déc 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')

Plus de réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by