Multiple plot on same graph

I would like to have these 4 plots on the same graph sheet for easy comparison. Plotting the individual plot works but I am not sure what is the reason why all the 4 plots cannot be plotted at the same time. Relying on benevolent members for assitance.

rng default
n = 30; 
% generate the A matrix
D = delsq(numgrid('L',n)); %same as matrix A
bGrid = ones(length(D),1); %same as RHS b vector
M = diag(diag(D));
tol = 1e-8;
maxit = 100;
P = inv(M);
L = ichol(M);
Q = ichol(M,struct('michol','on'));

% pcg  
[x0,f0,r0,it0,rv0] = pcg(D,bGrid); %no precond
[x1,f1,r1,it1,rv1] = pcg(D,bGrid,tol,maxit,P); % inv
[x2,f2,r2,it2,rv2] = pcg(D,bGrid,tol,maxit,L,L'); % ichol
[x3,f3,r3,it3,rv3] = pcg(D,bGrid,tol,maxit,Q,Q'); % modified ichol

% plots on same graph
semilogy(0:length(rv0)-1,rv0/norm(bGrid),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(bGrid),'-o')
semilogy(0:length(rv2)-1,rv2/norm(bGrid),'-o')
semilogy(0:length(rv3)-1,rv3/norm(bGrid),'-o')
yline(tol,'r--');
legend('No Preconditioner', 'inv(M)', 'ichol(M)', 'michol(M)', 'Tolerance', 'Location','southwest')
title('Conjugant Gredient Algorithm');
legend('boxoff')

Answer 1

[Walter Roberson]

rng default
n = 30; 
% generate the A matrix
D = delsq(numgrid('L',n)); %same as matrix A
bGrid = ones(length(D),1); %same as RHS b vector
M = diag(diag(D));
tol = 1e-8;
maxit = 100;
P = inv(M);
L = ichol(M);
Q = ichol(M,struct('michol','on'));

% pcg  
[x0,f0,r0,it0,rv0] = pcg(D,bGrid); %no precond
[x1,f1,r1,it1,rv1] = pcg(D,bGrid,tol,maxit,P); % inv
[x2,f2,r2,it2,rv2] = pcg(D,bGrid,tol,maxit,L,L'); % ichol
[x3,f3,r3,it3,rv3] = pcg(D,bGrid,tol,maxit,Q,Q'); % modified ichol

% plots on same graph
semilogy(0:length(rv0)-1,rv0/norm(bGrid),'-or')
hold on
semilogy(0:length(rv1)-1,rv1/norm(bGrid),'--.g')
semilogy(0:length(rv2)-1,rv2/norm(bGrid),':^b')
semilogy(0:length(rv3)-1,rv3/norm(bGrid),'+k')
yline(tol,'r--');
legend('No Preconditioner', 'inv(M)', 'ichol(M)', 'michol(M)', 'Tolerance', 'Location','southwest')
title('Conjugant Gredient Algorithm');
legend('boxoff')

min(rv0),max(rv0)

ans = 1.8506

ans = 52.6747

min(rv1),max(rv1)

ans = 1.9799e-07

ans = 52.6747

min(rv2),max(rv2)

ans = 1.9799e-07

ans = 52.6747

min(rv3),max(rv3)

ans = 1.9799e-07

ans = 52.6747

You are asking to solve the same system of linear equations with different preconditioners, but you are expecting that the differences between the results will be visible on the plot.

There is a notable difference on the first point of rv0 versus the remainder, but on the scale of the plot you cannot tell by the plot.

[Hmm!]

That is strange to know. What can be done about it? A similar plot here motivated my codes and all the plots were correctly diplayed. As indicated before, plotting individual plot yielded different shapes and positions and thus, combining all the 4 plots I expected to see all the plots not overlapping as it is in your case. Theoritically, I think there should be variations in the effectiveness of the different preconditioners.

[Walter Roberson]

They all interpolate to exactly the same thing, no difference at all.

rng default
n = 30; 
% generate the A matrix
D = delsq(numgrid('L',n)); %same as matrix A
bGrid = ones(length(D),1); %same as RHS b vector
M = diag(diag(D));
tol = 1e-8;
maxit = 100;
P = inv(M);
L = ichol(M);
Q = ichol(M,struct('michol','on'));

% pcg  
[x0,f0,r0,it0,rv0] = pcg(D,bGrid); %no precond
[x1,f1,r1,it1,rv1] = pcg(D,bGrid,tol,maxit,P); % inv
[x2,f2,r2,it2,rv2] = pcg(D,bGrid,tol,maxit,L,L'); % ichol
[x3,f3,r3,it3,rv3] = pcg(D,bGrid,tol,maxit,Q,Q'); % modified ichol

tx = [length(rv0), length(rv1), length(rv2), length(rv3)];
maxx = max(tx)

maxx = 68

xq = linspace(1, maxx, 1000);

whos

  Name         Size              Bytes  Class     Attributes

  D          588x588             49960  double    sparse    
  L          588x588             14120  double    sparse    
  M          588x588             14120  double    sparse    
  P          588x588             14120  double    sparse    
  Q          588x588             14120  double    sparse    
  bGrid      588x1                4704  double              
  f0           1x1                   8  double              
  f1           1x1                   8  double              
  f2           1x1                   8  double              
  f3           1x1                   8  double              
  it0          1x1                   8  double              
  it1          1x1                   8  double              
  it2          1x1                   8  double              
  it3          1x1                   8  double              
  maxit        1x1                   8  double              
  maxx         1x1                   8  double              
  n            1x1                   8  double              
  r0           1x1                   8  double              
  r1           1x1                   8  double              
  r2           1x1                   8  double              
  r3           1x1                   8  double              
  rv0         21x1                 168  double              
  rv1         68x1                 544  double              
  rv2         68x1                 544  double              
  rv3         68x1                 544  double              
  tol          1x1                   8  double              
  tx           1x4                  32  double              
  x0         588x1                4704  double              
  x1         588x1                4704  double              
  x2         588x1                4704  double              
  x3         588x1                4704  double              
  xq           1x1000             8000  double              

bgn = norm(bGrid)

bgn = 24.2487

r0q = interp1(1:length(rv0), rv0./bgn, xq);
r1q = interp1(1:length(rv1), rv1./bgn, xq);
r2q = interp1(1:length(rv2), rv2./bgn, xq);
r3q = interp1(1:length(rv3), rv3./bgn, xq);

[rv0(1), rv1(1), rv2(1), rv3(1)]

ans = 1×4

24.2487 24.2487 24.2487 24.2487

[r0q(1), r1q(1), r2q(1), r3q(1)]

ans = 1×4

1 1 1 1

subplot(4,1,1)
semilogy(xq, r0q)
title('rv0')
subplot(4,1,2)
plot(xq, r1q-r0q)
title('rv1 - rv0')
subplot(4,1,3)
plot(xq, r2q-r0q);
title('rv2 - rv0')
subplot(4,1,4)
plot(xq, r3q-r0q);
title('rv3 - rv0')

[Hmm!]

Is that the case for what is here? Scroll down to see the plot.

[Walter Roberson]

M = diag(diag(D));

That is a sparse diagonal matrix whose diagonal is the constant 4.

L = ichol(M);

That is a sparse diagonal matrix whose diagonal is the constant 2.

Your arrays are "boring" -- the hints you are providing through the additional parameters are adding essentially no usable new information, so the outputs are all the same (except being different sizes in some cases.)