Code is displaying graph but not plotting any points
조회 수: 2 (최근 30일)
이전 댓글 표시
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
채택된 답변
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Graph and Network Algorithms에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!