Algae growing. Concentration curve problem

2024 1 23
1 个回答
更新时间：2024 2 2
3 次查看（30 天）

0 个投票

Hi! i have made the following code:
function w = Bioreactor
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
alpha1 = 5*10^(-4); % Diffusion coefficient [m^2/s]
mu0 = 5;
mue = 3.5;
Hl = 5;
T = 5000; % [s]
nz = 100; nt = 1000;
KP = 2;
HP= 7;
PC = 1;
K = 2;
KN = 3;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5;
lambda = 2 ;
rs= 10;
P(1)=3; % Initial chosen concentration of P [g/L]
N(1)=9; % Initial chosen concentration of N [g/L]
C(1)=30; % Initial chosen concentration of C [g/L]
Z = 10; % Lenght of the bioreactor [m]
%Check of the problem's parameters
if any([P(1)-PC N(1)-NC] <=0 )
error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt); % Round nt to an integer
nz = round(nz); % Round nz to an integer
N = zeros(1, nt+1); % Initialize N to zero
P = zeros(1, nt+1); % Initialize P to zero
C = zeros(1, nt+1); % Initialize C to zero
I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);
% Finite-difference method
for j = 1:nt
f1= (KP*(P(j)-PC))/(HP+(P(j)-PC));
f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
pH = (6.35 - log10(C(j))) /2 ;
f3 = 1/(1+exp(lambda*(pH-PHs3)));
f4 = 1/(1+exp(lambda*(pH-PHs4)));
I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
g = mu0*I_in./(Hl+I_in);
integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j)); % Data organization
N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j)); % Function of N over the time and space
P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j)); % Function of P over the time and space
C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j)); % Function of C over the time and space
w(nz+1,j+1) = 9; % Algae concentration at z=0 [Dirichlet]
w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2);
w(1,j+1) = w(2,j+1); % Algae concentration at z=10 [Neumann]
plot(z,flipud(w(:,j+1)),'b');
%set(gca,'XDir','reverse');
xlabel('z'); ylabel('w');
title(['t=',num2str(t(j)),' s']);
axis([0 10 0 10]);
pause(0.0001);
end
end
And the code is working perfectly. But the W which is the algae concentration over space and time it has been fixed a value of 9 mg/L at the beginning of this bioreactor, so at z=0, im talking about the following row:
w(nz+1,j+1) = 9; % Algae concentration at z=0 [Dirichlet]
But now i would like to fix this concentration at a certain z which has to be equal to z=Z/2 so at the half of the bioreactor. I know that in order to make the code working for sure i have to modify this row:
w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2);
Can you help me to do that? I've tried something but anything really worked.
Thanks, attached to this post there is the file from which i took the equations.
I should have a graph like this:

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2024-1-23
编辑：Torsten 2024-1-23
Ran in:

0 个投票

Ran in:
Here is a code for a simple heat-conduction equation
dT/dt = d^2T/dz^2
T(z,0) = 0
T(0,t) = 200
T(1,t) = 340
T(0.5,t) = 90
You should be able to adapt it to your application.
L = 1;
zcut = L/2;
n1 = 50;
n2 = 50;
z1 = linspace(0,zcut,n1);
dz1 = z1(2)-z1(1);
z2 = linspace(zcut,L,n2);
dz2 = z2(2)-z2(1);
tstart = 0;
tend = 1;
dt = 1e-5;
nt = round((tend-tstart)/dt);
t = linspace(tstart,tend,nt);
T_at_zcut = 90;
T_at_0 = 200;
T_at_L = 340;
T = zeros(n1+n2-1,nt);
T(1,:) = T_at_0;
T(end,:) = T_at_L;
T(n1,:) = T_at_zcut;
for i = 1:nt-1
T(2:n1-1,i+1) = T(2:n1-1,i) + dt/dz1^2*(T(3:n1,i)-2*T(2:n1-1,i)+T(1:n1-2,i));
T(n1+1:n1+n2-2,i+1) = T(n1+1:n1+n2-2,i) + dt/dz2^2*(T(n1+2:n1+n2-1,i)-...
2*T(n1+1:n1+n2-2,i)+T(n1:n1+n2-3,i));
end
plot([z1,z2(2:end)],[T(:,1),T(:,100),T(:,250),T(:,500),T(:,750),T(:,1000),T(:,2500),T(:,5000),T(:,7500),T(:,end)])

26 个评论
显示 24更早的评论 隐藏 24更早的评论

Alfonso
Alfonso 2024-1-24
The following row:
T(n1,:) = T_at_zcut;
Is the one that fix the temperature in the center?
Torsten
Torsten 2024-1-24
编辑：Torsten 2024-1-24
Yes. The partial differential equation for T is only applied in i=2,...,n1-1,n1+1,...,n1+n2-2. The temperatures in the two boundary points and in the center are kept constant and are excluded from updates over time.
This is equivalent to solving two problems separately:
Problem 1:
dT/dt = d^2T/dz^2 0 <= z <= 0.5
T(z,0) = 0
T(0,t) = 200
T(0.5,t) = 90
Problem 2:
dT/dt = d^2T/dz^2 0.5 <= z <= 1
T(z,0) = 0
T(0.5,t) = 90
T(1,t) = 340
and glue the solutions at the right resp. left boundary points together.
Alfonso
Alfonso 2024-1-24
Ohh yes, i got this code. What do you think? I think it's working!
function Bioreactor_central
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
alpha1 = 0.005;
mu0 = 5;
mue = 3.5;
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2;
HP = 7;
PC = 1;
K = 2;
KN = 3;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5;
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
Z = 10;
% Check of the problem's parameters
if any([P(1)-PC N(1)-NC] <= 0 )
error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;
% Initialization
z = linspace(0, Z, nz+1);
t = linspace(0, T, nt+1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt+1);
P = zeros(1, nt+1);
C = zeros(1, nt+1);
I0 = @(t) 100*max(sin((t-0.4*3600)/(6.28)), sin(0));
w = zeros(nz+1, nt+1);
w(50,:) = 5;
% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1+n2-1, nt+1);
% Finite-difference method
for j = 1:nt
f1 = (KP*(P(j)-PC))/(HP+(P(j)-PC));
f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
pH = (6.35 - log10(C(j))) /2 ;
f3 = 1/(1+exp(lambda*(pH-PHs3)));
f4 = 1/(1+exp(lambda*(pH-PHs4)));
I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
g = mu0*I_in./(Hl+I_in);
integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));
N(j+1) = N(j) + dt*(-1/Z*integral_term*N(j));
P(j+1) = P(j) + dt*(-1/Z*integral_term*P(j));
C(j+1) = C(j) + dt*(-1/Z*integral_term*C(j));
w(nz+1,j+1) = w(nz,j+1);
w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + alpha1 * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + alpha1 * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
w(1,j+1) = w(2,j+1);
plot(z,w(:,j+1),'b');
xlabel('z'); ylabel('w');
title(['t=',num2str(t(j)),' s']);
axis([0 10 0 10]);
pause(0.0001);
end
end
Do you think that the results can be associated to the curve that i draw? I think so
Torsten
Torsten 2024-1-24
编辑：Torsten 2024-1-24
The setting
w(nz+1,j+1) = w(nz,j+1);
seems to be out of order because at the position of the code where this line appears, w(nz,j+1) is still set to 0 from the line of initialization
w = zeros(nz+1, nt+1);
So your boundary conditions for w are
dw/dx(0) = 0, w(L) = 0, w(L/2) = 5
Is this the correct setting you intend ?
Alfonso
Alfonso 2024-1-25
编辑：Alfonso 2024-1-25
i would like to have the same boundaries in the two border. I would like to have a neumann but not equal to 0. I would like to have the same behavior of the curve near the borders. I dont want 0.
I mean, in the code there is the light and it has to be put on a single side of the bioreactor, like the side at z=0. While on the other side, at z=10 i don't want the light so that with this plot i will see that from the central value, i have a bigger curve near the side with the light and a smaller curve near the side withtout the light.
Torsten
Torsten 2024-1-25
If you want dw/dz = 0 at both ends, use
w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + alpha1 * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + alpha1 * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);
instead of
w(nz+1,j+1) = w(nz,j+1);
w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + alpha1 * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + alpha1 * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
w(1,j+1) = w(2,j+1);
And since you decouple the two regions at z = L/2, I'm not sure if your integrations with trapz and cumtrapz still have to be over the complete region [0 L] or over [0 L/2] and [L/2 L] separately.
Alfonso
Alfonso 2024-1-25
编辑：Alfonso 2024-1-26
Ok @Torsten, i modified the code as you told me and so i get this code:
function Bioreactor_central
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5;
mue = 2;
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2;
HP = 7;
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5;
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
Z = 10;
% Check of the problem's parameters
if any([P(1) - PC N(1) - NC] <= 0 )
error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;
% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);
I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;
% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);
% Finite-difference method
for j = 1:nt
f1 = (KP * (P(j) - PC)) / (HP + (P(j) - PC));
f2 = (KN * (N(j) - NC)) / (HN + (N(j) - NC));
pH = (6.35 - log10(C(j))) / 2 ;
f3 = 1 / (1 + exp(lambda * (pH - PHs3)));
f4 = 1 / (1 + exp(lambda * (pH - PHs4)));
I_in = I0(t(j)) .* exp(-K * z.') .* exp((-rs) .* cumtrapz(z.', w(:, j)));
g = mu0 * I_in ./(Hl + I_in);
integral_term = f1 * f2 * f3 * trapz(z.', g .* w(:, j));
N(j + 1) = N(j) + dt * (-1 / Z * integral_term * N(j));
P(j + 1) = P(j) + dt * (-1 / Z * integral_term * P(j));
C(j + 1) = C(j) + dt * (-1 / Z * integral_term * C(j));
w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);
plot(z, w(:, j + 1), 'b');
xlabel('z'); ylabel('w');
title(['t=', num2str(t(j)), ' s']);
axis([0 10 0 10]);
pause(0.0001);
end
end
As a plot i have for example:
That seems correct because from the center value, i see 2 different curves with 2 different values on the borders. But, the I_in, which is the irradiation, is it put only on 1 side? Because at least i see that the 2 curves are increasing a lot, yes with 2 different values of w, but they are increasing together and i don't want something like this. In this case the plot should have a curve that has a permanent lower values of W and the other curve which is the oppost, i mean this other curve has to increase the w value with a different velocity due to the fact that the algae are irradiated by the light.
Torsten
Torsten 2024-1-26
I cannot help you to check whether your equations are adequate to model what you want to model. I can only help you to check whether predefined equations and boundary conditions are correctly implemented in MATLAB code.
Alfonso
Alfonso 2024-1-26
Ok, so. The equation of the irradiation is:
I_in = I0(t(j)) .* exp(-K * z.') .* exp((-rs) .* cumtrapz(z.', w(:, j)));
Now, i would like to know, if you can help me, where it is put. I mean i want the irradiation only on 1 side of the rectangule boundary, so it is like this? Or the irradiation is everywhere for all the z values?
Torsten
Torsten 2024-1-26
I cannot answer this.
At the moment, I_in is defined for all values of z and enters all equations for all values of z. Maybe you should save and plot the profiles to see what is happening.
Alfonso
Alfonso 2024-1-28
What do you mean that i should save and plot the profiles, which one, i didn't understand
Torsten
Torsten 2024-1-28
编辑：Torsten 2024-1-28
From what you write it seems to me that I_in as defined does not fit your needs.
So you could save I_in as you save w in a 2d-array and plot it to see if something goes wrong.
Alfonso
Alfonso 2024-1-31
@Torsten
I have the following updated code:
function Bioreactor_central
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5;
mue = 2;
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2;
HP = 7;
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5;
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
Z = 10;
% Check of the problem's parameters
if any([P(1) - PC N(1) - NC] <= 0 )
error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;
% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);
I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;
% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);
% Finite-difference method
for j = 1:nt
f1 = (KP * (P(j) - PC)) / (HP + (P(j) - PC));
f2 = (KN * (N(j) - NC)) / (HN + (N(j) - NC));
pH = (6.35 - log10(C(j))) / 2 ;
f3 = 1 / (1 + exp(lambda * (pH - PHs3)));
f4 = 1 / (1 + exp(lambda * (pH - PHs4)));
I_in = I0(t(j)) .* exp(-K * z.') .* exp((-rs) .* cumtrapz(z.', w(:, j)));
g = mu0 * I_in ./(Hl + I_in);
integral_term = f1 * f2 * f3 * trapz(z.', g .* w(:, j));
N(j + 1) = N(j) + dt * (-1 / Z * integral_term * N(j));
P(j + 1) = P(j) + dt * (-1 / Z * integral_term * P(j));
C(j + 1) = C(j) + dt * (-1 / Z * integral_term * C(j));
w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);
plot(z, w(:, j + 1), 'b');
xlabel('z'); ylabel('w');
title(['t=', num2str(t(j)), ' s']);
axis([0 10 0 10]);
pause(0.0001);
end
end
I understood that I_in is fixed for every value of the z axis, but instead i would like to have it only for like z=10 for every J value. I've tried to modify it and i got:
function Bioreactor_central
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5;
mue = 2;
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2;
HP = 7;
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5;
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
Z = 10;
% Check of the problem's parameters
if any([P(1) - PC N(1) - NC] <= 0 )
error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;
% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);
I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;
% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);
I_in = zeros(nz + 1, nt + 1);
for j = 1:nt
I_in(:, j) = I0(t(j)) * exp(-K * z.') .* exp((-rs) * cumtrapz(z, w(:, j).')');
end
% Finite-difference method
for j = 1:nt
f1 = (KP * (P(j) - PC)) / (HP + (P(j) - PC));
f2 = (KN * (N(j) - NC)) / (HN + (N(j) - NC));
pH = (6.35 - log10(C(j))) / 2 ;
f3 = 1 / (1 + exp(lambda * (pH - PHs3)));
f4 = 1 / (1 + exp(lambda * (pH - PHs4)));
I_in_10 = I_in(10, j);
g = mu0 * I_in_10 / (Hl + I_in_10);
integral_term = f1 * f2 * f3 * trapz(z, g .* w(:, j));
N(j + 1) = N(j) + dt * (-1 / Z * integral_term * N(j));
P(j + 1) = P(j) + dt * (-1 / Z * integral_term * P(j));
C(j + 1) = C(j) + dt * (-1 / Z * integral_term * C(j));
w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);
plot(z, w(:, j + 1), 'b');
xlabel('z'); ylabel('w');
title(['t=', num2str(t(j)), ' s']);
axis([0 10 0 10]);
pause(0.0001);
end
end
The problem is that the code works but doesn't give me anything as result. And i dont know why D:
Can you help me?
Torsten
Torsten 2024-1-31
Ran in:
Ran in:
Bioreactor_central()
function Bioreactor_central
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5 * 10^(-4);
mu0 = 5;
mue = 2;
Hl = 5;
T = 5000;
nz = 100; %nt = 1000;
KP = 2;
HP = 7;
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5;
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
Z = 10;
% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;
% Check of the problem's parameters
if any([P(1) - PC N(1) - NC] <= 0 )
error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
error('Check problem parameters')
end
% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);
I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;
I_in = zeros(nz + 1,nt);
% Finite-difference method
for j = 1:nt
f1 = (KP * (P(j) - PC)) / (HP + (P(j) - PC));
f2 = (KN * (N(j) - NC)) / (HN + (N(j) - NC));
pH = (6.35 - log10(C(j))) / 2 ;
f3 = 1 / (1 + exp(lambda * (pH - PHs3)));
f4 = 1 / (1 + exp(lambda * (pH - PHs4)));
I_in(:,j) = I0(t(j)) .* exp(-K * z.') .* exp((-rs) .* cumtrapz(z.', w(:, j)));
I_in_10 = I_in(Z,j);
g_in_10 = mu0 * I_in_10 / (Hl + I_in_10);
integral_term = f1 * f2 * f3 * trapz(z, g_in_10 .* w(:, j));
N(j + 1) = N(j) + dt * (-1 / Z * integral_term * N(j));
P(j + 1) = P(j) + dt * (-1 / Z * integral_term * P(j));
C(j + 1) = C(j) + dt * (-1 / Z * integral_term * C(j));
w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g_in_10 + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g_in_10 + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);
end
plot(z,w(:,end))
end
Alfonso
Alfonso 2024-2-1
编辑：Alfonso 2024-2-1
@Torsten
Ok, now the code is working again but i'm not having what i would like to have... :C
As you can see, i would like to have the red curve, not the blue which comes from the model. As i understood from the code that you sent the I_in is fixed only on z=10 and it is right but it doesnt happen what it should :/.
The algae should grow faster near the light side then the shadow side.
Alfonso
Alfonso 2024-2-1
@Torsten
I've tried to split everything and putting I_in equal to 0 from z=0 to z=5, so that i have the following code:
function Bioreactor_central
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5;
mue = 2;
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2;
HP = 7;
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5;
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
P_1(1) = 3;
N_1(1) = 9;
C_1(1) = 30;
Z = 10;
% Check of the problem's parameters
if any([P(1) - PC N(1) - NC] <= 0 )
error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;
% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);
N_1 = zeros(1, nt + 1);
P_1 = zeros(1, nt + 1);
C_1 = zeros(1, nt + 1);
I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;
% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);
% Finite-difference method
for j = 1:nt
f1 = (KP * (P(j) - PC)) / (HP + (P(j) - PC));
f2 = (KN * (N(j) - NC)) / (HN + (N(j) - NC));
pH = (6.35 - log10(C(j))) / 2 ;
f3 = 1 / (1 + exp(lambda * (pH - PHs3)));
f4 = 1 / (1 + exp(lambda * (pH - PHs4)));
f11 = (KP * (P_1(j) - PC)) / (HP + (P_1(j) - PC));
f22 = (KN * (N_1(j) - NC)) / (HN + (N_1(j) - NC));
pH1 = (6.35 - log10(C_1(j))) / 2 ;
f33 = 1 / (1 + exp(lambda * (pH1 - PHs3)));
f44 = 1 / (1 + exp(lambda * (pH1 - PHs4)));
I_in = zeros(nz + 1, 1);
I_in_1_5 = 0;%I0(t(j)) .* exp(-K * z(1:5).') .* exp((-rs) .* cumtrapz(z(1:5).', w(1:5, j)));
I_in_5_10 = I0(t(j)) .* exp(-K * z(5:10).') .* exp((-rs) .* cumtrapz(z(5:10).', w(5:10, j)));
g_1_5 = mu0 * I_in_1_5 ./(Hl + I_in_1_5);
g_5_10 = mu0 * I_in_5_10 ./(Hl + I_in_5_10);
integral_term_1_5 = f11 * f22 * f33 * trapz(z(1:5), g_1_5 .* w(1:5, j));
integral_term_5_10 = f1 * f2 * f3 * trapz(z(5:10), g_5_10 .* w(5:10, j));
N_1(j + 1) = N(j) + dt * (-1 / Z * integral_term_1_5 * N(j));
P_1(j + 1) = P(j) + dt * (-1 / Z * integral_term_1_5 * P(j));
C_1(j + 1) = C(j) + dt * (-1 / Z * integral_term_1_5 * C(j));
N(j + 1) = N(j) + dt * (-1 / Z * integral_term_5_10 * N(j));
P(j + 1) = P(j) + dt * (-1 / Z * integral_term_5_10 * P(j));
C(j + 1) = C(j) + dt * (-1 / Z * integral_term_5_10 * C(j));
w(2:49,j+1) = w(2:49,j) + dt * (f11*f22*f33*w(2:49,j).*g_1_5(2:49) + Dm * (w(3:50,j) - 2*w(2:49,j) + w(1:48,j))/dz^2) - mue*f4*w(2:49,j);
w(51:nz,j+1) = w(51:nz,j) + dt * (f1*f2*f3*w(51:nz,j).*g_5_10(51:nz) + Dm * (w(52:nz+1,j)-2*w(51:nz,j)+w(50:nz-1,j))/dz^2) - mue*f4*w(51:nz,j);
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);
plot(z, w(:, j + 1), 'b');
xlabel('z'); ylabel('w');
title(['t=', num2str(t(j)), ' s']);
axis([0 10 0 10]);
pause(0.0001);
end
%plot(z,w(end));
end
But i have the following error:
Index exceeds the number of array elements. Index must not exceed 1.
Error in Bioreactor_central1 (line 110)
w(2:49,j+1) = w(2:49,j) + dt * (f11*f22*f33*w(2:49,j).*g_1_5(2:49) + Dm * (w(3:50,j) - 2*w(2:49,j) + w(1:48,j))/dz^2) - mue*f4*w(2:49,j);
Do you know how can i solve?
Torsten
Torsten 2024-2-1
编辑：Torsten 2024-2-1
You split the domain and decouple all the equations over z=5. So why not making two computations on 0 <= z <= 5, once with I_in = 0 and second with the usual I_in ?
Alfonso
Alfonso 2024-2-1
编辑：Alfonso 2024-2-1
@Torsten
Yes, that's what i would like to do but it's hard for me to think about what to modify on matlab to obtain what i want xD
Torsten
Torsten 2024-2-1
Use your old code without the splitting of the w-array and the length of the region halfed - that's all.
Alfonso
Alfonso 2024-2-1
You are perfectly right! I've returned to the old code, the following:
function w = Bioreattore
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
alpha1 = 5*10^(-4);
mu0 = 5;
mue = 3.5;
Hl = 5;
T = 10000; % [s]
nz = 100; nt = 1000;
KP = 2 ;
HP= 7;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5;
lambda = 2 ;
rs= 10;
P(1)=3; % Initial chosen concentration of P
N(1)=9; % Initial chosen concentration of N
C(1)=30; % Initial chosen concentration of C
Z = 10; % Lenght of the bioreactor
%Check of the problem's parameters
if any([P(1)-PC N(1)-NC] <=0 )
error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt); % Round nt to an integer
nz = round(nz); % Round nz to an integer
N = zeros(1, nt+1); % Initialize N to zero
P = zeros(1, nt+1); % Initialize P to zero
C = zeros(1, nt+1); % Initialize C to zero
I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);
w(:,1) = 0; %Final concentration of the algae
% Finite-difference method
for j = 1:nt
f1= (KP*(P(j)-PC))/(HP+(P(j)-PC));
f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
pH = (6.35 - log10(C(j))) /2 ;
f3 = 1/(1+exp(lambda*(pH-PHs3)));
f4 = 1/(1+exp(lambda*(pH-PHs4)));
I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
g = mu0*I_in./(Hl+I_in);
integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j)); % Data organization
N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j)); % Function of N over the time and space
P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j)); % Function of P over the time and space
C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j)); % Function of C over the time and space
w(nz+1,j+1) = 5; % Algae concentration at z=0
w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2)- mue*f4*w(2:nz,j);
w(1,j+1) = w(2,j+1); % Algae concentration at z=10
plot(z,flipud(w(:,j+1)),'b');
%set(gca,'XDir','reverse');
xlabel('z'); ylabel('w');
title(['t=',num2str(t(j))]);
axis([0 10 0 10]);
pause(0.0001);
end
end
But at least while i remove put I_in equal to 0 i have the same plot, as follows:
What it's you advise?
Torsten
Torsten 2024-2-1
My advice remains the same:
Plot
I_in=I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
as you plot w to see if it's equal to or almost equal to 0.
Alfonso
Alfonso 2024-2-2
@Torsten
Ok, so. I did it, the code is the following:
function w = Bioreattore_Matlab
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% U(0,t) = 0
% U(L,t) = u (end-1,t)
% U(z, 0) = U0(z)=0
% it is assumed that the concentration of C,N,P are costant
% Problem parameters
alpha1 = 5*10^(-4);
mu0 = 5;
mue = 3.5;
Hl = 5;
T = 10000; % [s]
nz = 100; nt = 1000;
KP = 2 ;
HP= 7;
%P = 3;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
%N = 4;
NC = 2;
KC = 5;
HC = 5 ;
%C = 3 ;
PHs3 = 4 ;
PHs4 = 5;
Z = 10;
lambda = 2 ;
rs= 10;
P(1)=3;
N(1)=9;
C(1)=30;
%Check of the problem's parameters
if any([P(1)-PC N(1)-NC] <=0 )
error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T nz-2 nt-2 KP P PC HC KN N NC KC HC C HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz=Z/nz ;
%dt = T/nt;
%r = alpha1*dt/((dz)^2);
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;
%while (2*r) > st
% nt = nt+1;
%end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt); % Arrotonda nt a un numero intero
N = zeros(1, nt+1); % Inizializzazione di N a zero
P = zeros(1, nt+1); % Inizializzazione di P a zero
C = zeros(1, nt+1); % Inizializzazione di C a zero
%f1= (KP*(P-PC))/(HP+(P-PC));
%f2 = (KN*(N-NC))/(HN+(N-NC));
%pH = (6.35 - log10(C)) /2 ;
%f3 = 1/(1+exp(lambda*(pH-PHs3)));
%f4 = 1/(1+exp(lambda*(pH-PHs4)));
%I0 = 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
nz = round(nz);
nt = round(nt);
w = zeros(nz+1,nt+1);
w(:,1) = 0; %Final concentration of the algae [0]
% Finite-difference method
%for j=2:nt+1
% u(1) = 0; % Condizione al contorno di Dirichlet
% g = (mu0.*I(j))/(Hl+I(j)); decay = g(j)*f1*f2*f3-mue*f4;
% u(2:end-1)=(1-2*r-decay*dt)*u(2:end-1)+r*(u(1:end-2)+u(3:end));
% u(end) = u(end); % Condizione al contorno di Neumann
for j = 1:nt
f1= (KP*(P(j)-PC))/(HP+(P(j)-PC));
f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
pH = (6.35 - log10(C(j))) /2 ;
f3 = 1/(1+exp(lambda*(pH-PHs3)));
f4 = 1/(1+exp(lambda*(pH-PHs4)));
I_in = I0(t(j)).*exp(-K*z.').*exp((-rs));%.*cumtrapz(z.',w(:,j)));
g = mu0*I_in./(Hl+I_in);
integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));
N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j))
P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j))
C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j))
%I_in = I0(t(j))*exp(-K*z.').*exp(-rs)*cumtrapz(z.',w(:,j));
%w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz,j)+r*(w(1:nz-1,j)+w(3:nz+1,j));
w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2)
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = 5;%w(nz,j+1); %Concentration at z=0
%w(2:nz,j+1) = w(2:nz,j) + dt*...
%w(2:nz,j+1) = integral_term(j) +
plot(z,flipud(I_in),'b');
%set(gca,'XDir','reverse');
xlabel('z'); ylabel('w');
title(['t=',num2str(t(j))]);
%axis([0 10 0 10]);
pause(0.0001);
end
%plot(z,flipud(w(:,end)))
end
But the plots are not good as i expected, look at that:
They are the only plot that i see for all the time t is like 1 second the first plot and the next second the other. What do you think? What's you advise, they are just the I_in plot, dont read on yaxes the label.
Alfonso
Alfonso 2024-2-2
@Torsten
Regarding the old code:
function w = Bioreattore
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
alpha1 = 5*10^(-4);
mu0 = 0;
mue = 99999999;%3.5;
Hl = 5;
T = 15000; % [s]
nz = 100; nt = 1000;
KP = 2 ;
HP= 7;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5;
lambda = 2 ;
rs= 10;
P(1)=3; % Initial chosen concentration of P
N(1)=9; % Initial chosen concentration of N
C(1)=30; % Initial chosen concentration of C
Z = 10; % Lenght of the bioreactor
%Check of the problem's parameters
%if any([P(1)-PC N(1)-NC] <=0 )
% error('Check problem parameters')
%end
%if any([alpha1 mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
% error('Check problem parameters')
%end
% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt); % Round nt to an integer
nz = round(nz); % Round nz to an integer
N = zeros(1, nt+1); % Initialize N to zero
P = zeros(1, nt+1); % Initialize P to zero
C = zeros(1, nt+1); % Initialize C to zero
I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);
w(:,1) = 0; %Final concentration of the algae
% Finite-difference method
for j = 1:nt
f1= (KP*(P(j)-PC))/(HP+(P(j)-PC));
f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
pH = (6.35 - log10(C(j))) /2 ;
f3 = 1/(1+exp(lambda*(pH-PHs3)));
f4 = 1/(1+exp(lambda*(pH-PHs4)));
I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
g = mu0*I_in./(Hl+I_in);
integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j)); % Data organization
N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j)); % Function of N over the time and space
P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j)); % Function of P over the time and space
C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j)); % Function of C over the time and space
w(nz+1,j+1) = 5; % Algae concentration at z=0
w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2)- mue*f4*w(2:nz,j);
w(1,j+1) = w(2,j+1); % Algae concentration at z=10
plot(z,flipud(w(:,j+1)),'b');
%set(gca,'XDir','reverse');
xlabel('z'); ylabel('w');
title(['t=',num2str(t(j))]);
axis([0 10 0 10]);
pause(0.0001);
end
end
I see that if i change every parameter it doesn't change anything in the curve, itìs not normal, it seems that the w equation doesnt read the initial parameter, it gives me everytime the same curve also if i put everything 0. I just need a code that reads those changes at least, can you check?
Torsten
Torsten 2024-2-2
编辑：Torsten 2024-2-2
The parameters are all correctly used, but they are most probably far off. Did you take care about the units? In the article, some of them have time unit "day" and others have time unit "second". Did you convert them correctly to a common unit in your code ?
And the setting
w(nz+1,j+1) = 5;
is at the wrong position.
It should be set before the for-loop as
w(nz+1,:) = 5
Alfonso
Alfonso 2024-2-2
Ok, man i will check. At least i didn't check the dimension before D:
Alfonso
Alfonso 2024-2-5
@Torsten also with the right dimension the code doesn't work very well. Can i ask you a question? If i have this code:
function w = Bioreattore
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5*10^(-4); %m^2/s
mu0 = 8000%2.5; % L/day
mue = 3.5; % 1/day
Hl = 5; %half-saturation intensity W/(m^2*day)
T = 10000; % [s]
nz = 100; nt = 1000;
KP = 2 ;
HP= 7; %half-saturation concentrations of Phosphates
PC = 1; %critical concentration of the phosphates dopo la quale la crescita di w va a 0.
K = 2;
KN = 3 ;
HN = 0.02; %14.5*10^(-6); %half-saturation concentrations of nitrates [mgN/L]
NC = 2; %critical concentration of the nitrates dopo la quale la crescita di w va a 0.
KC = 5;
HC = 5 ;
PHs3 = 9 ; %describes the "switching" value of pH at which the growth increases if all other factors are kept unchanged.
lambda = 2 ; %sharpness of the profile of the f3 function
rs= 10; %L*m/d
P(1)=50; % Initial chosen concentration of P [mg/L]
N(1)=330; % Initial chosen concentration of N [mg/L]
C(1)=2100; % Initial chosen concentration of C/ Aumentando C, si abbassa il pH(+acido) e aumenta la morte algale [mg/L]
Z = 10; % Lenght of the bioreactor
%Check of the problem's parameters
if any([P(1)-PC N(1)-NC] <=0 )
error('Check problem parameters')
end
if any([Dm mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
error('Check problem parameters')
end
% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/Dm;
nt = T/dt;
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt); % Round nt to an integer
nz = round(nz); % Round nz to an integer
N = zeros(1, nt+1); % Initialize N to zero
P = zeros(1, nt+1); % Initialize P to zero
C = zeros(1, nt+1); % Initialize C to zero
I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);
w(:,1) = 0; %Final concentration of the algae
% Finite-difference method
for j = 1:nt
f1= (KP*(P(j)-PC))/(HP+(P(j)-PC))
f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC))
pH = (6.35 - log10(C(j))) /2 %relation between the pH and the concentration of carbon dioxide
f3 = 1/(1+exp(lambda*(pH-PHs3))) %The growth rate dependence function. Si riduce se il tutto diventa più acido (pH scende)..
Ha = mue*f3; %death rate of the algae. [gCOD/L]
I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j))); % W/m^2*s
g = mu0*I_in./(Hl+I_in);
integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j)); % Data organization
N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j)); % Function of N over the time and space
P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j)); % Function of P over the time and space
C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j)); % Function of C over the time and space
w(nz+1,j+1) = 5; % Algae concentration at z=0 [gSST/L]
w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + Dm * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2)- Ha*w(2:nz,j);
w(1,j+1) = w(2,j+1); % Algae concentration at z=10
plot(z,flipud(w(:,j+1)),'b');
%set(gca,'XDir','reverse');
xlabel('z'); ylabel('w');
title(['t=',num2str(t(j))]);
axis([0 10 0 10]);
pause(0.0001);
end
end
How can i modify the plot to obtain the following graphs:

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

标签

提问:

Alfonso
Alfonso
2024-1-23

评论：

Alfonso
Alfonso
2024-2-2

Community Treasure Hunt

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

Start Hunting!

Translated by