Итерационные вычисления в Matlab

Question

Я собираюсь смоделировать осаждение минерала во времени по ставке r. Когда этот минерал выпадает в осадок, он снижает концентрацию 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 против времени. Приведенный выше код не работает, потому что 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