I want to plot a numerical integral in a mesh which is dependent of x and y. How can I do this?

I have a two dimensional function and I want to plot the integral function in 2D mesh. I want to plot the C2 in color mesh and also want to plot diagonal element of C2. I am not getting the correct result. Please see the code I am working.

if true
  % 
L       =  20;
z1      =  0:0.01:L;
z2      =  0:0.01:L;
dz      = z1(2)-z1(1);
[Z1,Z2] = meshgrid(z1,z2);
W1    =   1+Z1/L;              
W2    =   2+Z2/L;              
for j   = 1:size(z1,2)
  for k= 1:size(z2,2)
      A(j,k) = sum(sum(1./W1(1:j,1:k))*dz)*dz;
      B(j,k) = sum(sum(1./W2(1:j,1:k))*dz)*dz;
  end
end
C2_h  =  (1-A).^2+(1-B);
figure(32);pcolor(Z1,Z2,C2_h);axis equal tight; shading flat;colorbar;colormap jet;
figure(5);
plot(z1/L, diag(C2_h) ,'r','Linewidth',2 );hold on;
end

Answer 1

[Anton Semechko]

The integrals in your equation can be evaluated analytically. Here how you can visualize C(z1,z2) and C(z1,z1):

L  = 20;
dz = 0.01;
z  = 0:dz:L;
[Z1,Z2] = meshgrid(z);

I1=@(z) log(1+z);
I2=@(z) log(2+z);
C=(1-I1(Z1)).^2 + (1+log(2)-I2(Z2));

figure('color','w')

ha=subplot(1,2,1);

z_lim=[z(1) z(end)];
imagesc(z_lim,z_lim,C)
axis equal 
hold on 

plot(z_lim(:),z_lim(:),'--k','LineWidth',2)
set(ha,'XLim',z_lim+[-1 1],'YLim',z_lim+[-1 1],'YDir','normal')
xlabel('z_1','Interpreter','tex','FontSize',20)
ylabel('z_2','Interpreter','tex','FontSize',20)

colorbar('Location','NorthOutside')

subplot(1,2,2)
Cd=diag(C);
plot(z(:),Cd(:),'--k','LineWidth',2)
xlabel('z_1','Interpreter','tex','FontSize',20)
ylabel('C(z_1,z_1)','Interpreter','tex','Interpreter','tex','FontSize',20)

[mk_ballav]

This is one of the special case where W1 and W2 are linear functions but I have other cases where W1 and W2 are not linear and I can't directly evaluate integral anlytically,so I have to do numerical integration. Here C(z1,z2) is the whole matrix elements and C(z1,z1) is just the diagonal elements of C.

[Anton Semechko]

OK. Assuming W1 and W2 are univariate functions of z1 and z2, respectively, you can integrate 1/W1 and 1/W2 numerically using builtin 'integral' function.. Here is an example:

% Visualization grid
L  = 20;
dz = 0.01;
z  = 0:dz:L;

% Evaluate definite integrals of 1/W1 and 1/W2 beween z(i-1) and z(i)
W1=@(z) 1./(z+1);
W2=@(z) 1./(z+2);
[I1,I2]=deal(zeros(size(z)));
for i=2:numel(z)
    I1(i)=integral(W1,z(i-1),z(i));
    I2(i)=integral(W2,z(i-1),z(i));
end

% Definite intergral between z(1) and z(i) = sum(I(1:i)) 
I1=cumsum(I1); 
I2=cumsum(I2);

% Function 
C=bsxfun(@plus,(1-I1).^2,(1-I2(:)));

% Visualization
figure('color','w')

ha=subplot(1,2,1);
z_lim=[z(1) z(end)];
imagesc(z_lim,z_lim,C)
axis equal
hold on

plot(z_lim(:),z_lim(:),'--k','LineWidth',2)
set(ha,'XLim',z_lim+[-1 1],'YLim',z_lim+[-1 1],'YDir','normal')
xlabel('z_1','Interpreter','tex','FontSize',20)
ylabel('z_2','Interpreter','tex','FontSize',20)

colorbar('Location','NorthOutside')

subplot(1,2,2)
Cd=diag(C);
plot(z(:),Cd(:),'--k','LineWidth',2)
xlabel('z_1','Interpreter','tex','FontSize',20)
ylabel('C(z_1,z_1)','Interpreter','tex','Interpreter','tex','FontSize',20)

You can confirm that the final result is visually indistinguishable from the one obtained analytically.