Variable layer thickness using loops

8 vues (au cours des 30 derniers jours)
Kathleen
Kathleen le 15 Sep 2023
Modifié(e) : Voss le 16 Sep 2023
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);

Réponse acceptée

Voss
Voss le 15 Sep 2023
% 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
M = 8
N = 3
N = 3
dz = [ 1 : 1 : N ] %variable layer thickness (m)
dz = 1×3
1 2 3
zt = cumsum( dz -dz(1) ) %layer tops (m)
zt = 1×3
0 1 3
zb = zt + dz %layer bottoms (m)
zb = 1×3
1 3 6
zmax = zb(end) %depth of model domain (m)
zmax = 6
rng(6)
sd = zmax * sort( rand(1,M) ) %sensor depths (m)
sd = 1×8
0.2502 0.6459 1.9919 2.5128 3.1789 3.5703 4.9274 5.3572
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 =========
%% 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
i = 1 jj = 1 i = 2 jj = 1 i = 3 jj = 2 i = 4 jj = 2 i = 5 jj = 3 i = 6 jj = 3 i = 7 jj = 3 i = 8 jj = 3
fprintf('\n======== single loop logical-vec method =========\n\n')
======== single loop logical-vec method =========
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 commentaires
Kathleen
Kathleen le 16 Sep 2023
Thank you so much!! I shed many tears over this assignment but I think I understand it now. Thank you so much!
Voss
Voss le 16 Sep 2023
Modifié(e) : Voss le 16 Sep 2023
You're welcome! Any questions, let me know.

Connectez-vous pour commenter.

Plus de réponses (1)

Walter Roberson
Walter Roberson le 15 Sep 2023
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

Catégories

En savoir plus sur Orange dans Help Center et File Exchange

Produits


Version

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by