Matlab中的迭代计算

Question

我想模拟一种矿物在时间上的沉淀速率，随着这种矿物的沉淀，它降低了溶液中Fe(II)的浓度，这降低了溶液相对于矿物的饱和状态，减慢了矿物沉淀的速度。

cla;
clf;

%%
%variables and rate calculation
Fe(1)=0.01*0.26/1000; %mM, * activity coefficient, convert to M
Si=1.8/1000; %mM, convert to M
pH=8;
H=10^(-pH);
k=1.11*10^(-9);
b=1.81;
Ksp=4.68*10^(22);
IAP=(Fe^(3)*Si^(2))/H^(6);
logIAP=log10(IAP);
logKsp=log10(Ksp); 

rate=k*exp(b*(logIAP-logKsp)); %mm Fe Kg-1 hour-1

%% 
%Fe vs time 
dt=0.01;
t=0:dt:10;

for ind=2:length(t)
Fe=Fe(1)-(k*exp(b*(logIAP-logKsp))/1000*t);
end

我最终想要一张Fe对time的曲线图。上面的代码不起作用，因为Fe依赖于logIAP，而后者又依赖于前一个时间步中的Fe ...我需要Matlab为每个时间步迭代地计算这两个值。

感谢您的帮助，谢谢！

Answer 1

正如其他人指出的那样，您需要使Fe成为一个数组。此外，您还需要在您的等式中正确索引t，因为这是一个向量。

Fe=zeros(size(t)); % now Fe is an array
Fe(1)=0.01*0.26/1000; % first value assigned
...
Ksp=4.68*10^(22);
IAP=(Fe(1)^(3)*Si^(2))/H^(6);
logIAP=log10(IAP);
logKsp=log10(Ksp); 
....
for ind=2:length(t)

    logIAP=log10((Fe(ind-1)^(3)*Si^(2))/H^(6));
    Fe(ind)=Fe(ind-1)-(k*exp(b*(logIAP-logKsp))/1000*t(ind));
end

Answer 2

使用Fe作为数组：

Fe=[0.01*0.26/1000]; % Fe is now an array

....

for ind=2:length(t)
    Fe(ind)=Fe(ind-1)-(k*exp(b*(logIAP-logKsp))/1000*t); % We fill Fe
end

然后将Fe绘制成t的函数