Is your code from above intended to solve equation (1) of the paper ?
Im working with algae growth linked to the light. I don't know how to write the integral on matlab inside the the light intensity formulation.
1 次查看（过去 30 天）
显示 更早的评论
function u = bioreactors
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 = 0.5;
mue = 0.035;
Hl = 5;
T = 10^3; % [s]
nz = 100; nt = 100;
KP = 0.7 ;
HP= 7;
P = 0.2 ;
PC = 0.1;
K = 0.50;
KN = 0.35 ;
HN = 14.5*10^(-6);
N = 0.3;
NC = 0.2;
KC = 0.5;
HC = 5 ;
C = 0.4 ;
PHs3 = 0.4 ;
PHs4 = 0.5;
Z = 10;
lambda = 2 ;
%Check of the problem's parameters
if any([P-PC N-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;
while (2*r) > st
nt = nt+1;
end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
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));
%I = ???% Dont know how to write it on matlab.
% 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
plot(u,z,'r*:');
set(gca,'YDir','reverse');
xlabel('u'); ylabel('z');
title(['t=',num2str(t(j))]);
axis([0 0.05 0 10]);
pause(0.01);
end
The integral is in the page 9 of the attachment "algae". The light intensity(I) formulation can be found at page 2, formulation number (3), it is a value which changes through the space z and the time.
As a result i should see a plot in which the mass of microorganism changes over the space and the time, something like this:
On Z=0 there is the light while in Z=5 there is no light so that's why the mass of microorganisms is equal to 0, they are photosynthetic algae.
回答（2 个）
Torsten
2023-11-18
移动：Torsten
2023-11-18
Where is your discretization of d^2w/dz^2 ? Where do you implement the boundary conditions dw/dz = 0 at z = 0 and z = L ?
You have to discretize equation (1) of the article in space and use an ode-integrator (usually ode15s) to integrate the values in the grid points over time.
Look up "method-of-lines" for more details.
I suggest you do this first without the integral term. After you succeeded, add the three ordinary differential equations for N,P and C and use "trapz" to compute the integral term.
1 个评论
Alfonso
2023-11-18
编辑：Alfonso
2023-11-18
i've already done the discretization of the equation (1) which is u(2:end-1)=(1-2*r-decay*dt)*u(2:end-1)+r*(u(1:end-2)+u(3:end)); but inside this formula there is the "decay" term which has "g" which depends by the "I" term which has an integral that I don't know how to insert inside matlab, this is the problem.
For the boundary condition i put:
- For z=0 --> u(1)=0, Dirichlet
- For z=L --> u(end)=u(end-1), Noemann equal to 0
You can find everything inside the for cycle. In fact, I just used a for cycle which is simpler to define the discretized w over the time, i'm a student so i'm not so good in doing this work.
But the problem that i don't know how to solve is that the I changes its values also during the space and my for cycle takes into account only the time and not the space.
Regarding N, P, C i want first to make the code going, then i'll go through their equations, so for now you see just some values of them without any formulas.
Torsten
2023-11-19
移动：Torsten
2023-11-19
My suggestion:
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w(:,1) = ... % w at time t = 0
for j = 1:nt
I_in = I0(t(j))*exp(-k*z.').*exp(-rs)*cumtrapz(z.',w(:,j));
g = mu0*I_in./(H_l+I_in);
integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));
w(2:nz,j+1) = w(2:nz,j) + dt*...
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);
end
27 个评论
Alfonso
2023-11-19
w(2:nz,j+1) = w(2:nz,j) + dt*...
what do you mean for dt*...? It that a function?
Torsten
2023-11-19
编辑：Torsten
2023-11-19
It means that you should insert the terms "rhs" of your differential equation
dw/dt = rhs(w,P,N,C,z)
in the grid points 2:nz.
This equation reads in discretized form
w(2:nz,j+1) = w(2:nz,j) + dt*rhs
You wrote you did this in the line
u(2:end-1)=(1-2*r-decay*dt)*u(2:end-1)+r*(u(1:end-2)+u(3:end));
I prefer not to overwrite the results from the previous time steps - thus my u is of dimension (nz+1,nt+1) while yours is of dimension nz+1.
Alfonso
2023-11-20
Im trying to follow your instruction. But i have an error. For the following code:
function u = bioreactors
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 = 0.5;
mue = 0.035;
Hl = 5;
T = 10^3; % [s]
nz = 100; nt = 100;
KP = 0.7 ;
HP= 7;
P = 0.2 ;
PC = 0.1;
K = 0.50;
KN = 0.35 ;
HN = 14.5*10^(-6);
N = 0.3;
NC = 0.2;
KC = 0.5;
HC = 5 ;
C = 0.4 ;
PHs3 = 0.4 ;
PHs4 = 0.5;
Z = 10;
lambda = 2 ;
rs= 10;
%Check of the problem's parameters
if any([P-PC N-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;
while (2*r) > st
nt = nt+1;
end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
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));
w(:,1) = 0.5
% 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
% w at time t = 0
for j = 1:nt
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));
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);
w(2:nz,j+1) = w(2:nz,j) + dt*...
plot(w,z,'r*:');
set(gca,'YDir','reverse');
xlabel('u'); ylabel('z');
title(['t=',num2str(t(j))]);
axis([0 0.05 0 10]);
pause(0.01);
end
I have the following error:
Dimension argument must be a positive integer scalar within indexing range.
Error in cumtrapz>getDimArg (line 91)
dim = matlab.internal.math.getdimarg(dim);
Error in cumtrapz (line 49)
dim = min(ndims(y)+1, getDimArg(dim));
Error in bioreactors (line 72)
I_in = I0(t(j))*exp(-K*z.').*exp(-rs)*cumtrapz(Z.',w(:,j));
Torsten
2023-11-20
w = zeros(nz+1,nt+1);
w(:,1) = 0.5
for j = 1:nt
I_in = I0(t(j)).*exp(-K*z.').*exp(-rs).*cumtrapz(z.',w(:,j));
...
Alfonso
2023-11-20
function u = bioreactors
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 = 0.5;
mue = 0.035;
Hl = 5;
T = 10^3; % [s]
nz = 100; nt = 100;
KP = 0.7 ;
HP= 7;
P = 0.2 ;
PC = 0.1;
K = 0.50;
KN = 0.35 ;
HN = 14.5*10^(-6);
N = 0.3;
NC = 0.2;
KC = 0.5;
HC = 5 ;
C = 0.4 ;
PHs3 = 0.4 ;
PHs4 = 0.5;
Z = 10;
lambda = 2 ;
rs= 10;
%Check of the problem's parameters
if any([P-PC N-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;
while (2*r) > st
nt = nt+1;
end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
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));
w = zeros(nz+1,nt+1);
w(:,1) = 0.5
% 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
I_in = I0(t(j)).*exp(-K*z.').*exp(-rs).*cumtrapz(z.',w(:,j));
%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));
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);
%w(2:nz,j+1) = w(2:nz,j) + dt*...
w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz,j+1)+r*(w(1:nz-1,j+1)+w(3:nz+1,j+1));
%w(2:nz,j+1) = integral_term(j) +
plot(w,z,'r*:');
set(gca,'YDir','reverse');
xlabel('w'); ylabel('z');
title(['t=',num2str(t(j))]);
axis([0 10 0 10]);
pause(0.01);
end
I tried to put the w(2:nz;j+1) in discretized mode but the graph doesnt work :C
Torsten
2023-11-20
w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz,j+1)+r*(w(1:nz-1,j+1)+w(3:nz+1,j+1));
You want to compute from t(j) to t(j+1). Thus the right-hand side can only have terms that are already computed. But you try to access w from time step j+1 that is not yet known.
And the plot command must be
plot(z,w(:,j+1))
instead of
plot(w,z,'r*:');
Alfonso
2023-11-20
Ok, i modified and the for cycle is now:
for j = 1:nt
I_in = I0(t(j)).*exp(-K*z.').*exp(-rs).*cumtrapz(z.',w(:,j));
%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));
w(1,j+1) = w(2,j+1);
w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz-1,j)+r*(w(1:nz-2,j)+w(3:nz,j));
w(nz+1,j+1) = w(nz,j+1);
Now it appears the problem:
Unable to perform assignment because the size of the left side is 99-by-1 and the size of the right
side is 98-by-1.
Error in bioreactors (line 78)
w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz-1,j)+r*(w(1:nz-2,j)+w(3:nz,j));
But i cannot find the problem
Torsten
2023-11-20
编辑：Torsten
2023-11-20
There are nz-1 elements from 2:nz. There are nz-2 elements from 3:nz and 2:nz-1 and nz-3 elements from 1:nz-2. Thus the assignment
w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz-1,j)+r*(w(1:nz-2,j)+w(3:nz,j));
is not possible.
My formula for I_in was not correct:
Replace
I_in = I0(t(j)).*exp(-K*z.').*exp(-rs).*cumtrapz(z.',w(:,j));
by
I_in = I0(t(j)).*exp(-K*z.').*exp(-rs*cumtrapz(z.',w(:,j));
And why do you continue to make wrong changes to my code suggestion ? The order of the assignments
w(2:nz,j+1) = w(2:nz,j) + dt*...
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);
matters.
w(1,j+1) = w(2,j+1);
w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz-1,j)+r*(w(1:nz-2,j)+w(3:nz,j));
w(nz+1,j+1) = w(nz,j+1);
makes no sense.
And by "rhs" I mean everything right from the = in
dt(w) = ...
(including the D_M*dzz(w) term).
Alfonso
2023-11-20
function u = bioreactors
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 = 10^4; % [s]
nz = 100; nt = 100;
KP = 2 ;
HP= 7;
P = 3 ;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
N = 3;
NC = 2;
KC = 5;
HC = 5 ;
C = 4 ;
PHs3 = 4 ;
PHs4 = 5;
Z = 10;
lambda = 2 ;
rs= 10;
%Check of the problem's parameters
if any([P-PC N-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;
while (2*r) > st
nt = nt+1;
end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
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));
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
I_in = I0(t(j)).*exp(-K*z.').*exp(-rs).*cumtrapz(z.',w(:,j));
%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));
w(1,j+1) = w(2,j+1);
%w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz-1,j)+r*(w(2:nz-3,j)+w(2:nz-2,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(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,w(:,j+1),'r*:');
%set(gca,'YDir','reverse');
xlabel('w'); ylabel('z');
title(['t=',num2str(t(j))]);
axis([0 10 0 10]);
pause(0.0001);
end
Thanks to you the code is working but there is a problem. I mean, while i put a value of T>1000 the code won't work, I mean, i don't see the plot. Which can be the problem?
Torsten
2023-11-20
bioreactors()
function u = bioreactors
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 = 5000; % [s]
nz = 100; nt = 1000;
KP = 2 ;
HP= 7;
P = 3 ;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
N = 3;
NC = 2;
KC = 5;
HC = 5 ;
C = 4 ;
PHs3 = 4 ;
PHs4 = 5;
Z = 10;
lambda = 2 ;
rs= 10;
%Check of the problem's parameters
if any([P-PC N-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;
%while (2*r) > st
%nt = nt+1;
%end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
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));
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
I_in = I0(t(j)).*exp(-K*z.').*exp(-rs*cumtrapz(z.',w(:,j)));
%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));
%w(1,j+1) = w(2,j+1);
%w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz-1,j)+r*(w(2:nz-3,j)+w(2:nz-2,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(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,w(:,j+1),'r*:');
%set(gca,'YDir','reverse');
%xlabel('w'); ylabel('z');
%title(['t=',num2str(t(j))]);
%axis([0 10 0 10]);
%pause(0.0001);
end
plot(z,w(:,end))
end
Alfonso
2023-11-20
编辑：Alfonso
2023-11-20
Ok, really thank you and sorry but im trying to do it for myself because my professor is not very good. You are awesome and the best for me! I've learned a lot talking with you then doing a 3 months course with my "prof"...
Regarding the plot, I would like that the plot starts from 0 on the axes and not from 10. I mean on a column 0 is supposed to be the top layer while z=10 is the bottom layer, so it's the one with no light and so without algae. But how can i change the code to do this?
Regarding the code that you sent to me, if you plot like this:
plot(z,w(:,j+1),'r*:');
xlabel('z'); ylabel('W');
title(['t=',num2str(t(j))]);
axis([0 10 0 10]);
pause(0.01);
You will see that the code won't work anymore. I mean, if i put like more then 5000 the curve stucks and won't change anymore during the time. Do you know why? Maybe, i have to check for the initial parameter?
Alfonso
2023-11-21
With this command it just reverse the axis but then the plot and the code starts always from 10. I want the code starting from 0, not from 10.
Torsten
2023-11-21
Could you sketch by hand what you want and include it here as a graphics ? I don't understand your description.
Alfonso
2023-11-21
编辑：Alfonso
2023-11-21
My code is the following:
function w = Bioreattore_POTENTE_(Mat)
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 = 1000; % [s]
nz = 100; nt = 100;
KP = 2 ;
HP= 7;
P = 3 ;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
N = 3;
NC = 2;
KC = 5;
HC = 5 ;
C = 4 ;
PHs3 = 4 ;
PHs4 = 5;
Z = 10;
lambda = 2 ;
rs= 10;
%Check of the problem's parameters
if any([P-PC N-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;
while (2*r) > st
nt = nt+1;
end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
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));
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
I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
%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));
%w(1,j+1) = w(2,j+1);
%w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz-1,j)+r*(w(2:nz-3,j)+w(2:nz-2,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(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,w(:,j+1),'r*:');
%set(gca,'XDir','reverse');
xlabel('z'); ylabel('w');
title(['t=',num2str(t(j))]);
axis([0 10 0 10]);
pause(0.0001);
end
%plot(z,w(:,end))
%end
Now, first of all my code works good until the time T=1000, after this value the code stops working and i don't know why. Second, the script starts at z=10 with a concentration value that i've chosen thanks to the following row in the for cycle:
w(1,j+1) = w(2,j+1); %Concentration at Z=0
w(nz+1,j+1) = 5; %Concentration at Z=10
But i want that the script starts at z=0 with a concentration value that i've chosen. The problem is that when i write the following row in the flow cycle:
w(1,j+1) = 5; %w(2,j+1); Concentration at Z=0
w(nz+1,j+1) = w(nz,j+1); %5; Concentration at Z=10
The script work but the plot not. I have the following, respectively, plots:
z=10; w=5
z=0; w=5Do you have some suggestions to solve my ploblems? Thank you very much
As attachment you can find the file that im following.
Torsten
2023-11-22
编辑：Torsten
2023-11-22
Your code works - also for T = 10000.
I already told you that the lines
%st = 1;
%while (2*r) > st
% nt = nt+1;
%end
don't make sense - so I left them, but as comments. They produce an endless loop if 2*r > st which is the case if T becomes large.
Further you use the Explicit Euler method for the discretization of d^2w/dz^2. The Explicit Euler method has strong stability problems of the time step size dt is too big. So if you get unphysical results, it's a good idea to choose a bigger value for nt.
Concerning your second question I don't know exactly what you want to do. Do you only want to reflect the graph at z = 5 or do you want to solve the equations with the boundary conditions reversed ? The latter will produce different results, I guess. If you only want to reflect results, you can use "flipud" on the results as done below.
w = Bioreattore_POTENTE_();
function w = Bioreattore_POTENTE_
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 = 2000;
KP = 2 ;
HP= 7;
P = 3 ;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
N = 3;
NC = 2;
KC = 5;
HC = 5 ;
C = 4 ;
PHs3 = 4 ;
PHs4 = 5;
Z = 10;
lambda = 2 ;
rs= 10;
%Check of the problem's parameters
if any([P-PC N-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;
%while (2*r) > st
% nt = nt+1;
%end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
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));
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
I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
%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));
%w(1,j+1) = w(2,j+1);
%w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz-1,j)+r*(w(2:nz-3,j)+w(2:nz-2,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(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,w(:,j+1),'r*:');
%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
Alfonso
2023-11-23
You delete everytime the plot that i put in the code, like the following:
%plot(z,w(:,j+1),'r*:');
%set(gca,'XDir','reverse');
%xlabel('z'); ylabel('w');
%title(['t=',num2str(t(j))]);
%axis([0 10 0 10]);
%pause(0.0001);
But i just want this type of plot, now the function flipud is good, thanks, but i would like to see how W changes over the time and the Z. If you try to run the code with the above plot you will see that after T=5000 the curve starts to freeze and stuck, the concentration of W doesn't change anymore after z=4 and it seems strange, the algae concentration should increase over the z and over the time but very slowly while we are far away from a light source, but as i told you in my above plot I see that the algae concentration near z=4 it stucks
Torsten
2023-11-23
You delete everytime the plot that i put in the code, like the following:
...
But i just want this type of plot, now the function flipud is good, thanks, but i would like to see how W changes over the time and the Z.
You can plot whatever you like after computing the solution up to time T. I would recommend separating calculations and postprocessing. But if you want to plot 100 or more curves within the for-loop: do it.
If you try to run the code with the above plot you will see that after T=5000 the curve starts to freeze and stuck, the concentration of W doesn't change anymore after z=4 and it seems strange, the algae concentration should increase over the z and over the time but very slowly while we are far away from a light source, but as i told you in my above plot I see that the algae concentration near z=4 it stucks
It indicates that your equation for w has reached a steady-state. If you think that no steady-state should exist, there must be something wrong with your equation. E.g. you did not yet consider P, N and C.
Alfonso
2024-1-24
编辑：Alfonso
2024-1-24
Ok, at least the code worked but now i noticed that regarding this code:
function w = Bioreactor
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5*10^(-4); % Diffusion coefficient [m^2/s]
mu0 = 5; % Bound for the growth rate
mue = 2;
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 = 0.2 ; %sharpness of the profile of f3
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([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);
% 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 [Dirichlet]
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)
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.1);
end
end
Which is the working one, but on the line:
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)
Was missing the parameter mue*f4 at the beginning of the code. So i added them like this:
w(2:nz,j+1) = mue*(f4*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)
But now if i click on run the code is instable, i tried to change nz and nt but without any solution, can you help me again @Torsten, you are the only one that can save me!
Torsten
2024-1-24
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 - mue*f4*w(2:nz,j))
Alfonso
2024-1-25
Thanks so much!. The only thing now is that the parameter mue*f4 is the death rate of the algae, so if i amplify mue there should be less algae and the curve should change a bit. But, trying to modify this death rate doesn't change the curve.
Can you check? Reducing and modifing the mue value?
Have you got some advises?
Torsten
2024-1-25
Have you got some advises?
No. The sink term is correctly implemented into the equation for w.
Alfonso
2024-1-25
Can you try to change the mue parameter? On the plot i see always the same curve :C
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)