I get an error that does not match the number of rows and columns in the program.

조회 수: 8 (최근 30일)
Pavel
Pavel 2023년 1월 20일
답변: Sulaymon Eshkabilov 2023년 1월 21일
I'm trying to execute code that I don't quite understand.
Error appears on startup
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in
the second matrix. To operate on each element of the matrix individually, use TIMES (.*) for elementwise multiplication.
Error in untitled>Xi2 (line 149)
M(:,(2*i+1):(2*i+2))=R*M(:,(2*(i-1)+1):(2*(i-1)+2));
Error in untitled (line 65)
Fmin=Xi2(zmin,ZJ,data,flux,kmed,W);
What I have to change? Can you explain or give links for study this problem?
function [] = untitled()
clear
flux=0.0167; % Atmospheric flux of 210Pb (Bq cm-2yr-1)
fulldepth=[1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35]; % mean depths of the layers of the sample (cm)
fulldata=[0.00844 0.00984 0.00777 0.01107 0.00921 0.01893 0.03159 0.03191 0.04810 0.03748 0.03277 0.02470 0.01529 0.01347 0.01016 0.00519 0.00197 0.00067]; % Mean 210Pb volumetric concentration in each layer of the sample (Bq cm-3)
zk=12; %Choice of depth of the washout
% We define new layers by gathering some of the upper ones just taking the corresponding averaged volumetric concentrations, listed in the vector 'data' below.
ZJ=[4,8,12,14,16,18,20,22,24,26,28,30,32,34,36]; %Vector of deepest points of each layer (i.e.divisions between layers)
data=[0.0091 0.0094 0.0141 0.03159 0.03191 0.04810 0.03748 0.03277 0.02470 0.01529 0.01347 0.01016 0.00519 0.00197 0.00067]; % mean 210Pb volumetric concentration profile gathering the upper sections (Bq cm-3)
W=ones(size(ZJ)); % Choice of the weighting of the fitting
kmax=size(ZJ,2);
i=1;
while(i<kmax)&&(ZJ(i)<zk)
i=i+1;
end
kmed=i;
j=kmed+1;
frac=data(j:size(data,2))/(trace(diag(data(j:size(data,2)))));
%% RANGE FOR THE OPTIMIZATION SEARCH
rmax=0.1; % Maximum value for the washout rates
Dmax=0.1; % Maximum value for the diffusive constant
wmax=0.48;% Maximum value for the acretion rate
wmin=0.24;% Minimum value for the acretion rate
%% INITIALIZATION OF THE SEARCH POINT
zmin=zeros(1,kmax+2);
zmin(j:kmax)=frac(1:(kmax-j+1));
zmin(1:kmed)=rmax*ones(1,kmed);
zmin(kmax+1)=wmax;
zmin(kmax+2)=Dmax;
z=zmin;
%% NUMBER OF INTERVALS IN EACH SEARCH DIRECTION
Nr=20; % number of intervals for the washout parameters
NW=10; % number of intervals for the advection parameter
ND=10; % number of intervals for the diffusive parameter
% Ordering of all the points (nr1, ..., nrkmed, nW, nD) with 1<= nri <= Nr, 1 <= nW <= NW, 1 <= nD<= ND, according to the assignment m=(nr1-1)*Nr^(kmed-1)*NW*ND + ... + (nrkmed-1)*NW*ND +(nW-1)*ND + nD, where 1<= m <= NW*ND*Nr^kmed
dim=Nr^kmed*NW*ND;
a=ceil((1:dim)./(NW*ND));
b=(1:dim)-(a-1)*NW*ND;
index=zeros(kmed+2,dim);
index(1,:)=ceil(a./Nr^(kmed-1));
for i=2:kmed
a=a-(index(i-1,:)-1).*Nr^(kmed-i+1);
index(i,:)=ceil(a./Nr^(kmed-i));
end
index(kmed+1,:)=ceil(b/ND);
index(kmed+2,:)=b-(index(kmed+1,:)-1)*ND;
clear a
clear b
disp(data)
disp(zmin)
disp(ZJ)
disp(flux)
disp(kmed)
disp(W)
Fmin=Xi2(zmin,ZJ,data,flux,kmed,W);
for i=1:dim
z(1:kmed)=rmax*index(1:kmed,i)/Nr;
z(kmax+1)=wmin+(wmax-wmin)*index(kmed+1,i)/NW;
z(kmax+2)=Dmax*index(kmed+2,i)/ND;
Faux=Xi2(z,ZJ,data,flux,kmed,W);
if(Faux<=Fmin)
Fmin=Faux;
zmin=z;
else
end
end
clear index
%% SAVING OF THE MATLAB FILE
save(['MFile_zk' num2str(ZJ(kmed)) '.mat'])
%% SAVING OF THE FIT RESULTS AND GRAPHICS PRINT OUT
f=figure;
Xmax=ZJ(size(ZJ,2))+4;
X=0:0.01:Xmax;
Ymax=0.2;
K=40;
C=Sol(zmin,X,ZJ,flux,kmed);
dx=X(2)-X(1);
imin=1;
imax=floor(Xmax/dx)+1;
step=70;
p=plot(X(imin:step:imax),C(imin:step:imax),'k',fulldepth,fulldata,'o','MarkerEdgeColor','k','MarkerFaceColor','k');
set(gca,'XTick',0:1:Xmax);
set(gca,'XTickLabel', ...
{'0','','','','','5','','','','','10','','','','','15','','','','','20','','','','','25','','','','','30','','','','','35','','','','','40'},'FontSize',15);
set(gca,'YTick',0:(Ymax/K):Ymax);
set(gca,'YTickLabel',{'0','','','','',num2str(5*Ymax/K),'','','','',num2str(10*Ymax/K),'','','','',num2str( ...
15*Ymax/K),'','','','',num2str(20*Ymax/K),'','','','',num2str(25*Ymax/K),'','','','',num2str(30*Ymax/K)},['' ...
'FontSize'],15);
saveas(p,['Fitting_zk' num2str(ZJ(kmed)) '.fig']);
close(f);
Inv=10^4*dx*sum((C(1:floor(ZJ(size(ZJ,2))/dx))+C(2:floor(ZJ(size(ZJ,2))/dx)+1))/2);
section_start=ZJ(1:kmed);
section_start(2:kmed)=ZJ(1:kmed-1);
section_start(1)=0;
fid=fopen(['Fitting_Result_zk' num2str(ZJ(kmed)) '.txt'],'w');
fprintf(fid,'\n%10s = %10.3f\n','flux (Bq cm-2 yr-1)',flux);
fprintf(fid,'\n%10s\n','ri (yr-1)');
fprintf(fid,'\n%10.3f over (%2.0f,%2.0f) cm\n',[zmin(1:kmed);section_start;ZJ(1:kmed)]);
fprintf(fid,'\n%10s\n\n','fi');
fprintf(fid,'%10.3f\n',frac);
fprintf(fid,'\n%10s = %1.3f\n','W(cm yr-1)',zmin(kmax+1));
fprintf(fid,'\n%10s = %1.3f\n','D(cm^2 yr-1)',zmin(kmax+2));
fprintf(fid,'\n%10s = %1.3f\n','X^2',Fmin);
fprintf(fid,'\n%10s = %1.3f\n','Inventory (Bq m-2)',Inv);
fprintf(fid,'\n%10s\n','Predicted Profile');
fprintf(fid,'\n%10s % 15s\n','depth(cm)','Pb210(Bq cm-3)');
fprintf(fid,'\n%10.1f % 15.4f',[X(imin:step:imax);C(imin:step:imax)]);
fclose(fid);
end
%% FUNCTIONS DEFINITIONS
function [f] = Xi2(X,Zk,Data,flux,kmed,W)
% GIVES THE VALUE OF THE SQUARED RESIDUALS FOR A GIVEN CONFIGURATION OF THE FITTING PARAMETERS
depth=Zk;
depthmed=Zk;
kmax=size(depth,2);
y=depth;
y(2:size(depth,2))=depth(2:size(depth,2))-depth(1:size(depth,2)-1);
lambda=0.0311;
thetaplus=zeros(1,kmed);
thetaminus=zeros(1,kmed);
Thetaplus=(X(kmax+1)+sqrt(X(kmax+1)^2+4*X(kmax+2)*lambda))/(2*X(kmax+2));
Thetaminus=(X(kmax+1)-sqrt(X(kmax+1)^2+4*X(kmax+2)*lambda))/(2*X(kmax+2));
M=zeros(2,2*kmax);
Ind=zeros(2,kmax);
for i=1:kmed
thetaplus(i)=(X(kmax+1)+sqrt(X(kmax+1)^2+4*X(kmax+2)*(lambda+X(i))))/(2*X(kmax+2));
thetaminus(i)=(X(kmax+1)-sqrt(X(kmax+1)^2+4*X(kmax+2)*(lambda+X(i))))/(2*X(kmax+2));
end
M(:,1:2)=[1 0;-thetaminus(1)/thetaplus(1) 0];
Ind(:,1)=[0;flux/(X(kmax+2)*thetaplus(1))];
for i=1:(kmed-1)
R=[exp(thetaplus(i)*depth(i))*(thetaplus(i)-thetaminus(i+1))/(exp(thetaplus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)))
exp(thetaminus(i)*depth(i))*(thetaminus(i)-thetaminus(i+1))/(exp(thetaplus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)));
exp(thetaplus(i)*depth(i))*(thetaplus(i+1)-thetaplus(i))/(exp(thetaminus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)))
exp(thetaminus(i)*depth(i))*(thetaplus(i+1)-thetaminus(i))/(exp(thetaminus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)))];
M(:,(2*i+1):(2*i+2))=R*M(:,(2*(i-1)+1):(2*(i-1)+2));
Ind(:,i+1)=R*Ind(:,i);
end
tS=X(1)*[(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
0 0]*M(:,1:2);
tT=X(1)*[(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
0 0]*Ind(:,1);
for j=2:kmed
tS=tS+X(j)*[(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j)
(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
0 0]*M(:,(2*(j-1)+1):(2*(j-1)+2));
tT=tT+X(j)*[(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j)
(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);0 0]*Ind(:,j);
end
M(:,(2*(kmax-1)+1):(2*(kmax-1)+2))=[0 0;0 1];
Ind(:,kmax)=[0;0];
for i=1:(kmax-(kmed+1))
S=X(1)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1)
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1)]*M(:,1:2);
T=X(1)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1)
(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1)
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1)]*Ind(:,1);
for j=2:kmed
S=S+X(j)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j)
(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j)
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j)]*M(:,(2*(j-1)+1):(2*(j-1)+2));
T=T+X(j)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j)
(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) (Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j)]*Ind(:,j);
end
M(:,(2*(kmax-i-1)+1):(2*(kmax-i-1)+2))=M(:,(2*(kmax-i)+1):(2*(kmax-i)+2))+(X(kmax-i)/y(kmax-i)-X(kmax-i+1)/y(kmax-i+1))*S/lambda;
Ind(:,kmax-i)=Ind(:,kmax-i+1)+(X(kmax-i)/y(kmax-i)-X(kmax-i+1)/y(kmax-i+1))*T/lambda;
end
tMkmed=[exp(thetaplus(kmed)*depth(kmed)) exp(thetaminus(kmed)*depth(kmed));
thetaplus(kmed)*exp(thetaplus(kmed)*depth(kmed))
thetaminus(kmed)*exp(thetaminus(kmed)*depth(kmed))];
tMkmed1=[exp(Thetaplus*depth(kmed)) exp(Thetaminus*depth(kmed));
Thetaplus*exp(Thetaplus*depth(kmed)) Thetaminus*exp(Thetaminus*depth(kmed))];
Mf=tMkmed*M(:,(2*(kmed-1)+1):(2*(kmed-1)+2))-tMkmed1*M(:,(2*kmed+1):(2*kmed+2))-(X(kmed+1)/y(kmed+1))*tS/lambda;
Nf=(X(kmed+1)/y(kmed+1))*tT/lambda+tMkmed1*Ind(:,kmed+1)-tMkmed*Ind(:,kmed);
MfInv=(1/(Mf(1,1)*Mf(2,2)-Mf(1,2)*Mf(2,1)))*[Mf(2,2) -Mf(1,2);-Mf(2,1) Mf(1,1)];
A=MfInv*Nf;
C=zeros(2,kmax);
for i=1:kmax
C(:,i)=M(:,(2*(i-1)+1):(2*(i-1)+2))*A+Ind(:,i);
end
F=0;
G=Data*diag(W)*transpose(Data);
for i=1:size(depthmed,2)
if (i<=kmed)
F = F+W(i)*(C(1,i)*exp(depthmed(i)*thetaplus(i))+C(2,i)*exp(depthmed(i)*thetaminus(i))-Data(i))^2;
else
B=(X(i)/y(i))*(tS*A+tT)/lambda;
if(C(1,i)==0)
F =F+ W(i)*(C(2,i)*exp(depthmed(i)*Thetaminus)+B(1)-Data(i))^2;
else
F =F+ W(i)*(C(1,i)*exp(depthmed(i)*Thetaplus)+C(2,i)*exp(depthmed(i)*Thetaminus)+B(1)-Data(i))^2;
end
end
end
f=F/G;
end
function [V] = Sol(X,Z,depth,flux,kmed)
% SOLVES THE BOUNDARY AND CONTINUITY CONDITIONS FOR A GIVEN
CONFIGURATION OF THE PARAMETERS
kmax=size(depth,2);
y=depth;
y(2:kmax)=depth(2:kmax)-depth(1:kmax-1);
lambda=0.0311;
thetaplus=zeros(1,kmed);
thetaminus=zeros(1,kmed);
Thetaplus=(X(kmax+1)+sqrt(X(kmax+1)^2+4*X(kmax+2)*lambda))/(2*X(kmax+2));
Thetaminus=(X(kmax+1)-sqrt(X(kmax+1)^2+4*X(kmax+2)*lambda))/(2*X(kmax+2));
M=zeros(2,2*kmax);
Ind=zeros(2,kmax);
for i=1:kmed
thetaplus(i)=(X(kmax+1)+sqrt(X(kmax+1)^2+4*X(kmax+2)*(lambda+X(i))))/(2*X(kmax+2));
thetaminus(i)=(X(kmax+1)-sqrt(X(kmax+1)^2+4*X(kmax+2)*(lambda+X(i))))/(2*X(kmax+2));
end
M(:,1:2)=[1 0;-thetaminus(1)/thetaplus(1) 0];
Ind(:,1)=[0;flux/(X(kmax+2)*thetaplus(1))];
for i=1:(kmed-1)
R=[exp(thetaplus(i)*depth(i))*(thetaplus(i)-thetaminus(i+1))/(exp(thetaplus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)))
exp(thetaminus(i)*depth(i))*(thetaminus(i)-thetaminus(i+1))/(exp(thetaplus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)));
exp(thetaplus(i)*depth(i))*(thetaplus(i+1)-thetaplus(i))/(exp(thetaminus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)))
exp(thetaminus(i)*depth(i))*(thetaplus(i+1)-thetaminus(i))/(exp(thetaminus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)))];
M(:,(2*i+1):(2*i+2))=R*M(:,(2*(i-1)+1):(2*(i-1)+2));
Ind(:,i+1)=R*Ind(:,i);
end
tS=X(1)*[(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
0 0]*M(:,1:2);
tT=X(1)*[(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
0 0]*Ind(:,1);
for j=2:kmed
tS=tS+X(j)*[(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j)
(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
0 0]*M(:,(2*(j-1)+1):(2*(j-1)+2));
tT=tT+X(j)*[(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j)
(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
0 0]*Ind(:,j);
end
M(:,(2*(kmax-1)+1):(2*(kmax-1)+2))=[0 0;0 1];
Ind(:,kmax)=[0;0];
for i=1:(kmax-(kmed+1))
S=X(1)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1)
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1)]*M(:,1:2);
T=X(1)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1)
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1)]*Ind(:,1);
for j=2:kmed
S=S+X(j)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j)
(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) (Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j)]*M(:,(2*(j-1)+1):(2*(j-1)+2));
T=T+X(j)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j)
(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) (Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j)]*Ind(:,j);
end
M(:,(2*(kmax-i-1)+1):(2*(kmax-i-1)+2))=M(:,(2*(kmax-i)+1):(2*(kmax-i)+2))+(X(kmax-i)/y(kmax-i)-X(kmax-i+1)/y(kmax-i+1))*S/lambda;
Ind(:,kmax-i)=Ind(:,kmax-i+1)+(X(kmax-i)/y(kmax-i)-X(kmax-i+1)/y(kmax-i+1))*T/lambda;
end
tMkmed=[exp(thetaplus(kmed)*depth(kmed)) exp(thetaminus(kmed)*depth(kmed));
thetaplus(kmed)*exp(thetaplus(kmed)*depth(kmed))
14
thetaminus(kmed)*exp(thetaminus(kmed)*depth(kmed))];
tMkmed1=[exp(Thetaplus*depth(kmed)) exp(Thetaminus*depth(kmed));
Thetaplus*exp(Thetaplus*depth(kmed)) Thetaminus*exp(Thetaminus*depth(kmed))];
Mf=tMkmed*M(:,(2*(kmed-1)+1):(2*(kmed-1)+2))-tMkmed1*M(:,(2*kmed+1):(2*kmed+2))-(X(kmed+1)/y(kmed+1))*tS/lambda;
Nf=(X(kmed+1)/y(kmed+1))*tT/lambda+tMkmed1*Ind(:,kmed+1)-tMkmed*Ind(:,kmed);
MfInv=(1/(Mf(1,1)*Mf(2,2)-Mf(1,2)*Mf(2,1)))*[Mf(2,2) -Mf(1,2);-Mf(2,1) Mf(1,1)];
A=MfInv*Nf;
C=zeros(2,kmax);
for i=1:kmax
C(:,i)=M(:,(2*(i-1)+1):(2*(i-1)+2))*A+Ind(:,i);
end
ConVal=zeros(1,size(Z,2));
j=1;
for i=1:size(Z,2)
while(j<kmax)&&(Z(i)>depth(j))
j=j+1;
end
if (j<=kmed)
ConVal(i)=C(1,j)*exp(Z(i)*thetaplus(j))+C(2,j)*exp(Z(i)*thetaminus(j));
else
B=(X(j)/y(j))*(tS*A+tT)/lambda;
if(C(1,j)==0)
ConVal(i)=C(2,j)*exp(Z(i)*Thetaminus)+B(1);
else
ConVal(i)=C(1,j)*exp(Z(i)*Thetaplus)+C(2,j)*exp(Z(i)*Thetaminus)+B(1);
end
end
end
V=ConVal;
end

이 질문에 답변하려면 로그인하십시오.

채택된 답변

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2023년 1월 21일
There were over dozens of size related errors, figure handle and uncommented comments in your code that are fixed. The main function is also renamed as RUNALL (untitled is not a good name :) for such a long code). Now it runs and completes the whole simulation.
function [] = RUNALL()
clear
flux=0.0167; % Atmospheric flux of 210Pb (Bq cm-2yr-1)
fulldepth=[1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35]; % mean depths of the layers of the sample (cm)
fulldata=[0.00844 0.00984 0.00777 0.01107 0.00921 0.01893 0.03159 0.03191 0.04810 0.03748 0.03277 0.02470 0.01529 0.01347 0.01016 0.00519 0.00197 0.00067]; % Mean 210Pb volumetric concentration in each layer of the sample (Bq cm-3)
zk=12; %Choice of depth of the washout
% We define new layers by gathering some of the upper ones just taking the corresponding averaged volumetric concentrations, listed in the vector 'data' below.
ZJ=[4,8,12,14,16,18,20,22,24,26,28,30,32,34,36]; %Vector of deepest points of each layer (i.e.divisions between layers)
data=[0.0091 0.0094 0.0141 0.03159 0.03191 0.04810 0.03748 0.03277 0.02470 0.01529 0.01347 0.01016 0.00519 0.00197 0.00067]; % mean 210Pb volumetric concentration profile gathering the upper sections (Bq cm-3)
W=ones(size(ZJ)); % Choice of the weighting of the fitting
kmax=size(ZJ,2);
i=1;
while(i<kmax)&&(ZJ(i)<zk)
i=i+1;
end
kmed=i;
j=kmed+1;
frac=data(j:size(data,2))/(trace(diag(data(j:size(data,2)))));
%% RANGE FOR THE OPTIMIZATION SEARCH
rmax=0.1; % Maximum value for the washout rates
Dmax=0.1; % Maximum value for the diffusive constant
wmax=0.48;% Maximum value for the acretion rate
wmin=0.24;% Minimum value for the acretion rate
%% INITIALIZATION OF THE SEARCH POINT
zmin=zeros(1,kmax+2);
zmin(j:kmax)=frac(1:(kmax-j+1));
zmin(1:kmed)=rmax*ones(1,kmed);
zmin(kmax+1)=wmax;
zmin(kmax+2)=Dmax;
z=zmin;
%% NUMBER OF INTERVALS IN EACH SEARCH DIRECTION
Nr=20; % number of intervals for the washout parameters
NW=10; % number of intervals for the advection parameter
ND=10; % number of intervals for the diffusive parameter
% Ordering of all the points (nr1, ..., nrkmed, nW, nD) with 1<= nri <= Nr, 1 <= nW <= NW, 1 <= nD<= ND, according to the assignment m=(nr1-1)*Nr^(kmed-1)*NW*ND + ... + (nrkmed-1)*NW*ND +(nW-1)*ND + nD, where 1<= m <= NW*ND*Nr^kmed
dim=Nr^kmed*NW*ND;
a=ceil((1:dim)./(NW*ND));
b=(1:dim)-(a-1)*NW*ND;
index=zeros(kmed+2,dim);
index(1,:)=ceil(a./Nr^(kmed-1));
for i=2:kmed
a=a-(index(i-1,:)-1).*Nr^(kmed-i+1);
index(i,:)=ceil(a./Nr^(kmed-i));
end
index(kmed+1,:)=ceil(b/ND);
index(kmed+2,:)=b-(index(kmed+1,:)-1)*ND;
clear a
clear b
disp(data)
disp(zmin)
disp(ZJ)
disp(flux)
disp(kmed)
disp(W)
Fmin=Xi2(zmin,ZJ,data,flux,kmed,W);
for i=1:dim
z(1:kmed)=rmax*index(1:kmed,i)/Nr;
z(kmax+1)=wmin+(wmax-wmin)*index(kmed+1,i)/NW;
z(kmax+2)=Dmax*index(kmed+2,i)/ND;
Faux=Xi2(z,ZJ,data,flux,kmed,W);
if(Faux<=Fmin)
Fmin=Faux;
zmin=z;
else
end
end
clear index
%% SAVING OF THE MATLAB FILE
save(['MFile_zk' num2str(ZJ(kmed)) '.mat'])
%% SAVING OF THE FIT RESULTS AND GRAPHICS PRINT OUT
f=figure;
Xmax=ZJ(size(ZJ,2))+4;
X=0:0.01:Xmax;
Ymax=0.2;
K=40;
C=Sol(zmin,X,ZJ,flux,kmed);
dx=X(2)-X(1);
imin=1;
imax=floor(Xmax/dx)+1;
step=70;
plot(X(imin:step:imax),C(imin:step:imax),'k',fulldepth,fulldata,'o','MarkerEdgeColor','k','MarkerFaceColor','k');
set(gca,'XTick',0:1:Xmax);
set(gca,'XTickLabel', ...
{'0','','','','','5','','','','','10','','','','','15','','','','','20','','','','','25','','','','','30','','','','','35','','','','','40'},'FontSize',15);
set(gca,'YTick',0:(Ymax/K):Ymax);
set(gca,'YTickLabel',{'0','','','','',num2str(5*Ymax/K),'','','','',num2str(10*Ymax/K),'','','','',num2str( ...
15*Ymax/K),'','','','',num2str(20*Ymax/K),'','','','',num2str(25*Ymax/K),'','','','',num2str(30*Ymax/K)},['' ...
'FontSize'],15);
saveas(f,['Fitting_zk' num2str(ZJ(kmed)) '.fig']);
close(f);
Inv=10^4*dx*sum((C(1:floor(ZJ(size(ZJ,2))/dx))+C(2:floor(ZJ(size(ZJ,2))/dx)+1))/2);
section_start=ZJ(1:kmed);
section_start(2:kmed)=ZJ(1:kmed-1);
section_start(1)=0;
fid=fopen(['Fitting_Result_zk' num2str(ZJ(kmed)) '.txt'],'w');
fprintf(fid,'\n%10s = %10.3f\n','flux (Bq cm-2 yr-1)',flux);
fprintf(fid,'\n%10s\n','ri (yr-1)');
fprintf(fid,'\n%10.3f over (%2.0f,%2.0f) cm\n',[zmin(1:kmed);section_start;ZJ(1:kmed)]);
fprintf(fid,'\n%10s\n\n','fi');
fprintf(fid,'%10.3f\n',frac);
fprintf(fid,'\n%10s = %1.3f\n','W(cm yr-1)',zmin(kmax+1));
fprintf(fid,'\n%10s = %1.3f\n','D(cm^2 yr-1)',zmin(kmax+2));
fprintf(fid,'\n%10s = %1.3f\n','X^2',Fmin);
fprintf(fid,'\n%10s = %1.3f\n','Inventory (Bq m-2)',Inv);
fprintf(fid,'\n%10s\n','Predicted Profile');
fprintf(fid,'\n%10s % 15s\n','depth(cm)','Pb210(Bq cm-3)');
fprintf(fid,'\n%10.1f % 15.4f',[X(imin:step:imax);C(imin:step:imax)]);
fclose(fid);
end
% FUNCTIONS DEFINITIONS
function [f] = Xi2(X,Zk,Data,flux,kmed,W)
% GIVES THE VALUE OF THE SQUARED RESIDUALS FOR A GIVEN CONFIGURATION OF THE FITTING PARAMETERS
depth=Zk;
depthmed=Zk;
kmax=size(depth,2);
y=depth;
y(2:size(depth,2))=depth(2:size(depth,2))-depth(1:size(depth,2)-1);
lambda=0.0311;
thetaplus=zeros(1,kmed);
thetaminus=zeros(1,kmed);
Thetaplus=(X(kmax+1)+sqrt(X(kmax+1)^2+4*X(kmax+2)*lambda))/(2*X(kmax+2));
Thetaminus=(X(kmax+1)-sqrt(X(kmax+1)^2+4*X(kmax+2)*lambda))/(2*X(kmax+2));
M=zeros(2,2*kmax);
Ind=zeros(2,kmax);
for i=1:kmed
thetaplus(i)=(X(kmax+1)+sqrt(X(kmax+1)^2+4*X(kmax+2)*(lambda+X(i))))/(2*X(kmax+2));
thetaminus(i)=(X(kmax+1)-sqrt(X(kmax+1)^2+4*X(kmax+2)*(lambda+X(i))))/(2*X(kmax+2));
end
M(:,1:2)=[1 0;-thetaminus(1)/thetaplus(1) 0];
Ind(:,1)=[0;flux/(X(kmax+2)*thetaplus(1))];
for i=1:(kmed-1)
R=[exp(thetaplus(i)*depth(i))*(thetaplus(i)-thetaminus(i+1))/(exp(thetaplus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1))) ...
exp(thetaminus(i)*depth(i))*(thetaminus(i)-thetaminus(i+1))/(exp(thetaplus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)));
exp(thetaplus(i)*depth(i))*(thetaplus(i+1)-thetaplus(i))/(exp(thetaminus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1))) ...
exp(thetaminus(i)*depth(i))*(thetaplus(i+1)-thetaminus(i))/(exp(thetaminus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)))];
M(:,(2*i+1):(2*i+2))=R*M(:,(2*(i-1)+1):(2*(i-1)+2));
Ind(:,i+1)=R*Ind(:,i);
end
tS=X(1)*[(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
0 0]*M(:,1:2);
tT=X(1)*[(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
0 0]*Ind(:,1);
for j=2:kmed
tS=tS+X(j)*[(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) ...
(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
0 0]*M(:,(2*(j-1)+1):(2*(j-1)+2));
tT=tT+X(j)*[(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) ...
(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);0 0]*Ind(:,j);
end
M(:,(2*(kmax-1)+1):(2*(kmax-1)+2))=[0 0;0 1];
Ind(:,kmax)=[0;0];
for i=1:(kmax-(kmed+1))
S=X(1)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1), (Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) ...
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1)]*M(:,1:2);
T=X(1)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) ...
(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1)...
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1)]*Ind(:,1);
for j=2:kmed
S=S+X(j)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) ...
(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) ...
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j)]*M(:,(2*(j-1)+1):(2*(j-1)+2));
T=T+X(j)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) ...
(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) (Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j)]*Ind(:,j);
end
M(:,(2*(kmax-i-1)+1):(2*(kmax-i-1)+2))=M(:,(2*(kmax-i)+1):(2*(kmax-i)+2))+(X(kmax-i)/y(kmax-i)-X(kmax-i+1)/y(kmax-i+1))*S/lambda;
Ind(:,kmax-i)=Ind(:,kmax-i+1)+(X(kmax-i)/y(kmax-i)-X(kmax-i+1)/y(kmax-i+1))*T/lambda;
end
tMkmed=[exp(thetaplus(kmed)*depth(kmed)), exp(thetaminus(kmed)*depth(kmed));
thetaplus(kmed)*exp(thetaplus(kmed)*depth(kmed)) ...
thetaminus(kmed)*exp(thetaminus(kmed)*depth(kmed))];
tMkmed1=[exp(Thetaplus*depth(kmed)) exp(Thetaminus*depth(kmed));
Thetaplus*exp(Thetaplus*depth(kmed)) Thetaminus*exp(Thetaminus*depth(kmed))];
Mf=tMkmed*M(:,(2*(kmed-1)+1):(2*(kmed-1)+2))-tMkmed1*M(:,(2*kmed+1):(2*kmed+2))-(X(kmed+1)/y(kmed+1))*tS/lambda;
Nf=(X(kmed+1)/y(kmed+1))*tT/lambda+tMkmed1*Ind(:,kmed+1)-tMkmed*Ind(:,kmed);
MfInv=(1/(Mf(1,1)*Mf(2,2)-Mf(1,2)*Mf(2,1)))*[Mf(2,2) -Mf(1,2);-Mf(2,1) Mf(1,1)];
A=MfInv*Nf;
C=zeros(2,kmax);
for i=1:kmax
C(:,i)=M(:,(2*(i-1)+1):(2*(i-1)+2))*A+Ind(:,i);
end
F=0;
G=Data*diag(W)*transpose(Data);
for i=1:size(depthmed,2)
if (i<=kmed)
F = F+W(i)*(C(1,i)*exp(depthmed(i)*thetaplus(i))+C(2,i)*exp(depthmed(i)*thetaminus(i))-Data(i))^2;
else
B=(X(i)/y(i))*(tS*A+tT)/lambda;
if(C(1,i)==0)
F =F+ W(i)*(C(2,i)*exp(depthmed(i)*Thetaminus)+B(1)-Data(i))^2;
else
F =F+ W(i)*(C(1,i)*exp(depthmed(i)*Thetaplus)+C(2,i)*exp(depthmed(i)*Thetaminus)+B(1)-Data(i))^2;
end
end
end
f=F/G;
end
function [V] = Sol(X,Z,depth,flux,kmed)
% SOLVES THE BOUNDARY AND CONTINUITY CONDITIONS FOR A GIVEN
% CONFIGURATION OF THE PARAMETERS
kmax=size(depth,2);
y=depth;
y(2:kmax)=depth(2:kmax)-depth(1:kmax-1);
lambda=0.0311;
thetaplus=zeros(1,kmed);
thetaminus=zeros(1,kmed);
Thetaplus=(X(kmax+1)+sqrt(X(kmax+1)^2+4*X(kmax+2)*lambda))/(2*X(kmax+2));
Thetaminus=(X(kmax+1)-sqrt(X(kmax+1)^2+4*X(kmax+2)*lambda))/(2*X(kmax+2));
M=zeros(2,2*kmax);
Ind=zeros(2,kmax);
for i=1:kmed
thetaplus(i)=(X(kmax+1)+sqrt(X(kmax+1)^2+4*X(kmax+2)*(lambda+X(i))))/(2*X(kmax+2));
thetaminus(i)=(X(kmax+1)-sqrt(X(kmax+1)^2+4*X(kmax+2)*(lambda+X(i))))/(2*X(kmax+2));
end
M(:,1:2)=[1 0;-thetaminus(1)/thetaplus(1) 0];
Ind(:,1)=[0;flux/(X(kmax+2)*thetaplus(1))];
for i=1:(kmed-1)
R=[exp(thetaplus(i)*depth(i))*(thetaplus(i)-thetaminus(i+1))/(exp(thetaplus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1))) ...
exp(thetaminus(i)*depth(i))*(thetaminus(i)-thetaminus(i+1))/(exp(thetaplus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)));
exp(thetaplus(i)*depth(i))*(thetaplus(i+1)-thetaplus(i))/(exp(thetaminus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1))) ...
exp(thetaminus(i)*depth(i))*(thetaplus(i+1)-thetaminus(i))/(exp(thetaminus(i+1)*depth(i))*(thetaplus(i+1)-thetaminus(i+1)))];
M(:,(2*i+1):(2*i+2))=R*M(:,(2*(i-1)+1):(2*(i-1)+2));
Ind(:,i+1)=R*Ind(:,i);
end
tS=X(1)*[(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
0 0]*M(:,1:2);
tT=X(1)*[(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
0 0]*Ind(:,1);
for j=2:kmed
tS=tS+X(j)*[(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) ...
(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
0 0]*M(:,(2*(j-1)+1):(2*(j-1)+2));
tT=tT+X(j)*[(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) ...
(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
0 0]*Ind(:,j);
end
M(:,(2*(kmax-1)+1):(2*(kmax-1)+2))=[0 0;0 1];
Ind(:,kmax)=[0;0];
for i=1:(kmax-(kmed+1))
S=X(1)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) (Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) ...
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1)]*M(:,1:2);
T=X(1)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1), (Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(1)*depth(1))-1)/thetaplus(1) ...
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(1)*depth(1))-1)/thetaminus(1)]*Ind(:,1);
for j=2:kmed
S=S+X(j)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) ...
(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j), (Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j)]*M(:,(2*(j-1)+1):(2*(j-1)+2));
T=T+X(j)*[(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j) ...
(Thetaminus*exp(-Thetaplus*depth(kmax-i))/(Thetaplus-Thetaminus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j);
(Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaplus(j)*depth(j))-exp(thetaplus(j)*depth(j-1)))/thetaplus(j), (Thetaplus*exp(-Thetaminus*depth(kmax-i))/(Thetaminus-Thetaplus))*(exp(thetaminus(j)*depth(j))-exp(thetaminus(j)*depth(j-1)))/thetaminus(j)]*Ind(:,j);
end
M(:,(2*(kmax-i-1)+1):(2*(kmax-i-1)+2))=M(:,(2*(kmax-i)+1):(2*(kmax-i)+2))+(X(kmax-i)/y(kmax-i)-X(kmax-i+1)/y(kmax-i+1))*S/lambda;
Ind(:,kmax-i)=Ind(:,kmax-i+1)+(X(kmax-i)/y(kmax-i)-X(kmax-i+1)/y(kmax-i+1))*T/lambda;
end
tMkmed=[exp(thetaplus(kmed)*depth(kmed)) exp(thetaminus(kmed)*depth(kmed));
thetaplus(kmed)*exp(thetaplus(kmed)*depth(kmed)) ...
thetaminus(kmed)*exp(thetaminus(kmed)*depth(kmed))];
tMkmed1=[exp(Thetaplus*depth(kmed)) exp(Thetaminus*depth(kmed));
Thetaplus*exp(Thetaplus*depth(kmed)) Thetaminus*exp(Thetaminus*depth(kmed))];
Mf=tMkmed*M(:,(2*(kmed-1)+1):(2*(kmed-1)+2))-tMkmed1*M(:,(2*kmed+1):(2*kmed+2))-(X(kmed+1)/y(kmed+1))*tS/lambda;
Nf=(X(kmed+1)/y(kmed+1))*tT/lambda+tMkmed1*Ind(:,kmed+1)-tMkmed*Ind(:,kmed);
MfInv=(1/(Mf(1,1)*Mf(2,2)-Mf(1,2)*Mf(2,1)))*[Mf(2,2) -Mf(1,2);-Mf(2,1) Mf(1,1)];
A=MfInv*Nf;
C=zeros(2,kmax);
for i=1:kmax
C(:,i)=M(:,(2*(i-1)+1):(2*(i-1)+2))*A+Ind(:,i);
end
ConVal=zeros(1,size(Z,2));
j=1;
for i=1:size(Z,2)
while(j<kmax)&&(Z(i)>depth(j))
j=j+1;
end
if (j<=kmed)
ConVal(i)=C(1,j)*exp(Z(i)*thetaplus(j))+C(2,j)*exp(Z(i)*thetaminus(j));
else
B=(X(j)/y(j))*(tS*A+tT)/lambda;
if(C(1,j)==0)
ConVal(i)=C(2,j)*exp(Z(i)*Thetaminus)+B(1);
else
ConVal(i)=C(1,j)*exp(Z(i)*Thetaplus)+C(2,j)*exp(Z(i)*Thetaminus)+B(1);
end
end
end
V=ConVal;
end

추가 답변 (0개)

카테고리

Help Center 및 File Exchange에서 Programmatic Model Editing에 대해 자세히 알아보기

태그

제품


릴리스

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by