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

8 次查看(过去 30 天)
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

回答(1 个)

Anton Semechko
Anton Semechko 2018-7-4
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)
  2 个评论
mk_ballav
mk_ballav 2018-7-4
编辑:mk_ballav 2018-7-4
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
Anton Semechko 2018-7-4
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.

请先登录,再进行评论。

Community Treasure Hunt

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

Start Hunting!

Translated by