Variable layer thickness using loops
이전 댓글 표시
Good afternoon,
I am trying to make the graphs as attached. However I keep running into error code:
Error using /
Matrix dimensions must agree.
Error in HW4_inverse3_KVW (line 46)
jj = ceil( sd / dz); %jj is layer that sd(i) resides within
dz is variable layer thickness dz = 1 2 3 (m).
How can I make a code for the variable layer thickness? I do not understand what is going on.
Thank you
% Make VSP G-mat for arbitrary sensor depths and non-uniform layer thicknesses
clear; close all; clc
set(0,'DefaultAxesLineWidth',2,'DefaultLineLineWidth',3)
set(0,'DefaultAxesXminorTick','on','DefaultAxesYminorTick','on')
set(0,'DefaultAxesFontSize',12,'DefaultAxesFontWeight','bold')
M = 8
N = 3
dz = [ 1 : 1 : N ] %variable layer thickness (m)
zt = cumsum( dz -dz(1) ) %layer tops (m)
zb = zt + dz %layer bottoms (m)
zmax = zb(end) %depth of model domain (m)
rng(6)
sd = zmax * sort( rand(1,M) ) %sensor depths (m)
figure(1)
plot(ones(1,M), sd,'bsq')
yline( [zt zb(end)], 'r-', 'linewidth',2 )
ylim([0 zmax]); set(gca,'ydir','reverse','Xtick', [ ] )
title('sensors and layers'); ylabel('depth (m)')
title('Sensors and Layers');
ylabel('Depth (m)')
G0 = zeros( M, N );
fprintf('======== single loop simplest method =========\n\n')
%% single loop simplest method (least lines)
for i = 1: M
jj = ceil( sd(i) / dz); %jj is layer that sd(i) resides within
G0( i, jj ) = sd(i) - zt( jj ); %assign segment length in jj layer
G0( i, 1: jj - 1) = dz; %works for jj=1 because 1:0 is not executed
fprintf('i = %d jj = %d \n', i, jj )
end
fprintf('\n======== single loop logical-vec method =========\n\n')
figure(2)
subplot (2,2,1)
imagesc(G0); colorbar; colormap('cool(20)'); axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G0')
G1 = zeros(M,N);
채택된 답변
% Make VSP G-mat for arbitrary sensor depths and non-uniform layer thicknesses
clear; close all; clc
set(0,'DefaultAxesLineWidth',2,'DefaultLineLineWidth',3)
set(0,'DefaultAxesXminorTick','on','DefaultAxesYminorTick','on')
set(0,'DefaultAxesFontSize',12,'DefaultAxesFontWeight','bold')
M = 8
N = 3
dz = [ 1 : 1 : N ] %variable layer thickness (m)
zt = cumsum( dz -dz(1) ) %layer tops (m)
zb = zt + dz %layer bottoms (m)
zmax = zb(end) %depth of model domain (m)
rng(6)
sd = zmax * sort( rand(1,M) ) %sensor depths (m)
figure(1)
plot(ones(1,M), sd,'bsq')
yline( [zt zb(end)], 'r-', 'linewidth',2 )
ylim([0 zmax]); set(gca,'ydir','reverse','Xtick', [ ] )
title('sensors and layers'); ylabel('depth (m)')
title('Sensors and Layers');
ylabel('Depth (m)')
G0 = zeros( M, N );
fprintf('======== single loop simplest method =========\n\n')
%% single loop simplest method (least lines)
for i = 1: M
% jj = ceil( sd(i) / dz ) %jj is layer that sd(i) resides within
jj = find(zt <= sd(i), 1, 'last');
G0( i, jj ) = sd(i) - zt( jj ); %assign segment length in jj layer
% G0( i, 1: jj - 1) = dz; %works for jj=1 because 1:0 is not executed
G0( i, 1: jj - 1) = dz( 1: jj - 1);
fprintf('i = %d jj = %d \n', i, jj )
end
fprintf('\n======== single loop logical-vec method =========\n\n')
figure(2)
% subplot (2,2,1)
imagesc(G0); colorbar; colormap('cool(20)'); axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G0')
댓글 수: 8
Walter Roberson
2023년 9월 15일
Right at the very top of Voss's posting there are some tool icons to link the posting or flag it, and so on.
On the line immediately below that in Voss's posting, to the right, it says "Ran in R2023b"
On the right hand side of the line below the "Ran in" there is a gray box followed by the word "Copy". Click that word "Copy". Now go over to your MATLAB editor and use your system Paste facility (for example control-V for Windows) to get a copy of Voss's code in the editor on your system.
Kathleen
2023년 9월 15일
Voss
2023년 9월 15일
The fundamental problem, I think, is that these code snippets treat dz as a scalar, when in fact it is a vector.
Kathleen
2023년 9월 15일
See below for changes required. The method for calculating G1 required indexing into dz rather than multiplying dz by a vector of ones (and assigning elements 1 through nl of row i on the left-hand side), and the method for calculating G2 just required replacing dz by dz(j).
% Make VSP G-mat for arbitrary sensor depths and non-uniform layer thicknesses
clear; close all; clc
set(0,'DefaultAxesLineWidth',2,'DefaultLineLineWidth',3)
set(0,'DefaultAxesXminorTick','on','DefaultAxesYminorTick','on')
set(0,'DefaultAxesFontSize',12,'DefaultAxesFontWeight','bold')
M = 8
N = 3
dz = [ 1 : 1 : N ] %variable layer thickness (m)
zt = cumsum( dz -dz(1) ) %layer tops (m)
zb = zt + dz %layer bottoms (m)
zmax = zb(end) %depth of model domain (m)
rng(6)
sd = zmax * sort( rand(1,M) ) %sensor depths (m)
figure(1)
plot(ones(1,M), sd,'bsq')
yline( [zt zb(end)], 'r-', 'linewidth',2 )
ylim([0 zmax]); set(gca,'ydir','reverse','Xtick', [ ] )
title('sensors and layers'); ylabel('depth (m)')
title('Sensors and Layers');
ylabel('Depth (m)')
%% single loop simplest method (fewest lines)
fprintf('======== single loop simplest method =========\n\n')
G0 = zeros( M, N );
for i = 1: M
% jj = ceil( sd(i) / dz ) %jj is layer that sd(i) resides within
jj = find(zt <= sd(i), 1, 'last');
G0( i, jj ) = sd(i) - zt( jj ); %assign segment length in jj layer
% G0( i, 1: jj - 1) = dz; %works for jj=1 because 1:0 is not executed
G0( i, 1: jj - 1) = dz( 1: jj - 1);
fprintf('i = %d jj = %d \n', i, jj )
end
figure(2)
subplot (2,2,1)
imagesc(G0); colorbar; colormap('cool(20)'); axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G0')
%% vectorized single loop logical-vec method
fprintf('\n======== single loop logical-vec method =========\n\n')
G1 = zeros(M,N);
for i = 1 : M
ip = sd(i) - zt >= 0; % logical: layers crossed by i'th ray == true (1)
nl = sum( ip ); % num layers crossed by ith ray
blrl = sd(i) - zt(nl); % bottom layer ray length
% G1( i ,ip ) = [ ones(1,nl-1)*dz blrl ];
G1( i , 1:nl ) = [ dz(1:nl-1) blrl ];
fprintf('i= %3d nl= %3d sd = %.4f blrl= %.4f\n', i,nl,sd(i), blrl )
end
subplot(2,2,2)
imagesc(G1); colorbar; colormap('cool(20)'); axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G1')
%% double loop method
fprintf('\n===========double loop method ======\n\n')
G2 = zeros(M,N);
for i = 1 : M
for j = 1 : N
d = sd( i ) - zt (j ); %distance between i'th sensor and j'th layer top
if d <= dz(j) & d >= 0 %sensor within j'th layer
G2( i , j ) = d;
elseif sd( i ) > zt( j ) %sensor crosses j'th layer
G2( i , j ) = dz(j);
end %if
fprintf('i=%3d j=%3d sd=%.2f zt=%.2f d= %+.2f G(i,j)=%.2f\n', ...
i,j,sd(i),zt(j), d, G2(i,j) )
end %inner loop
fprintf('\n')
end %outer loop
subplot(2,2,3)
imagesc(G2); colorbar; colormap('cool(20)'); axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G2')
Kathleen
2023년 9월 16일
추가 답변 (1개)
Walter Roberson
2023년 9월 15일
sd(i) / dz
sd(i) is a scalar but dz is a row vector.
In MATLAB, the / operator is mrdivide, / matrix right divide. A/B is similar to A * pinv(B) where * is algebraic matrix multiplication. For the / operator to work properly, the number of columns in the numerator and denominator must be the same . You have a scalar numerator, so one column, and your denominiator is a row vector with N columns. Unless N is 1, then the / operator would fail.
If you wanted [sd(i) / dz(1), sd(i) / dz(2), sd(i) / dz(3)] and so on, a collection of individual divisions, then you need to use rdivide, ./
sd(i) ./ dz
카테고리
도움말 센터 및 File Exchange에서 Loops and Conditional Statements에 대해 자세히 알아보기
Translated by 
