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.

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.

[Torsten]

Is your code from above intended to solve equation (1) of the paper ?

[Alfonso]

Yes, i want just to have the growth of microorganisms over the space and time. If you dont have some parameters just choose a value random, then i'll change it with real ones.

Answer 1

[Torsten]

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.

[Alfonso]

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.

Answer 2

[Torsten]

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

[Alfonso]

w(2:nz,j+1) = w(2:nz,j) + dt*...
    
what do you mean for dt*...? It that a function?

[Walter Roberson]

dt = t(2)-t(1)

[Alfonso]

I knew that, but i dont know what the "..." means. It gives me error on matlab, should i put something instead of ...?

[Torsten]

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]

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]

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]

ok, now i have to add the rhs terms. But are the following?

Thank you very much for the support and sorry for disturbing you a lot!

[Alfonso]

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]

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]

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]

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]

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]

You still have the two wrong settings in your code that I described in my last comment.

And the code doesn't work because 2*r > st, thus MATLAB never leaves the while-loop

while (2*r)  > st
nt = nt+1;
end

[Torsten]

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]

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?

[Walter Roberson]

@Alfonso perhaps you are looking for

set(gca, 'Xdir', 'reverse')

[Alfonso]

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]

Could you sketch by hand what you want and include it here as a graphics ? I don't understand your description.

[Alfonso]

@Torsten

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=5

Do you have some suggestions to solve my ploblems? Thank you very much

As attachment you can find the file that im following.

[Torsten]

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]

Hi @Torsten, thank you everytime for helping me.

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]

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]

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]

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]

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]

Have you got some advises?

No. The sink term is correctly implemented into the equation for w.

[Alfonso]

Can you try to change the mue parameter? On the plot i see always the same curve :C