Algae growing, concentration curve problem

Hi, i have made the following code:

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 = 8000; % [s]
nz = 100; nt = 2000;
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;
%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))); %Equation (3) in Algae
    %I_in = I0(t(j))*exp(-K*z.').*exp(-rs)*cumtrapz(z.',w(:,j));
    g = mu0*I_in./(Hl+I_in); %Equation (2) in Algae
    integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j)); %The first part of the equation (1)
    %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)); %Equation (1) in Algae
    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(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

It regards the growing of algae in a tower in which there is light only on 1 side of it. So, first of all, at z=0 we put some seeds so that we have a certain concentration of algae and then when i start the code i would like to see the path of the algae concentration (w) over the whole tower and over the time.

At this point, the code that i've made works really good until T=2000 (seconds) but after the code stops working and so the curve that i see on the printed graph is like "stucked" and so it wont work anymore.

Im loosing my mind to find some problems, but i really dont know what to modify, i've tried to do everything i know. Can you help me to solve si problem, i would like to reach also T=20000 seconds without any problem to my curve printed.

I really appriciate you hand, so please help me :C

P.S. Attached to this post you find also the pdf file of the article that i read to make the code.

[Sam Chak]

Hi @Alfonso

I noticed that you used cumtrapz() and trapz() in the code, and I'm checking if you can actually model the differential equations like that. Could you please edit and add comments, such as the equation number, in the code? This will help to verify if the equations can be effectively coded in MATLAB.

Additionally, it seems that the integral terms in Eq. (9) can be approximated by the sum, for i=1 to m=100.

Excerpt: Kumar, K., Pisarenco, M., Rudnaya, M., & Savcenco, V. (2014). A note on analysis and numerics of algae growth. Nonlinear Analysis : Real World Applications, 15(1), 392-403.

[Alfonso]

I've modified the code with some comments like Equation (3) in Algae.

Regarding the equation of f1, f2, f3, they are funcion, respectivetely of P, N and C. Now first i want that the code works, so i hypotized P, N and C as constants so alfo f1, f2 and f3 are constant. But while we solve the problem that i explained above, we should go through the implementation of the equation of P, N and C which changes over time and space like you can see in the equations (9) in the algae pdf.

Answer 1

[Alan Stevens]

Increase nt to 10000 and it works just fine for T = 20000.

[Alfonso]

i've done it but if you do this and wait ti.l the end at a certain point you see that the curve goes up and down but its stucked...

[Image Analyst]

Show the plot for us.

[Alan Stevens]

@Alfonso "i've done it but if you do this and wait ti.l the end at a certain point you see that the curve goes up and down but its stucked..." Really?!! This is what I get with what is basically your program:

% Problem parameters
alpha1 =  5*10^(-4); 
mu0 = 5; 
mue = 3.5; 
Hl = 5;
T = 20000; % [s]
nz = 100; nt = 10000;
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;
%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;

% If the following condition were ever invoked it would loop forever!!
% 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]
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(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
end
plot(z,flipud(w(:,end))),grid
xlabel('z'); ylabel('w');
title(['t=',num2str(t(end))]); 
axis([0 10 0 10]);

[Sam Chak]

Hi @Alfonso

Are you expecting these plots?

[Alfonso]

Nono, im not expecting them, maybe in the future ahahahha. What firstly i need is the graph that you obtained, but for example if you fix T=40000 you will see that the curve stucks at a certain point, so for example also at T=40000 you will have the same graph of T=20000. I mean if you see the concentration of the algaes (W) at z=6 it will be always 0, also if i fix T equal to 50000 or 100000 and this is strange.

[Alfonso]

https://www.dropbox.com/scl/fi/hitrucftdf7xams238kri/Plot_problem.mp4?rlkey=oxxyhqhlgclt3vyrsgr8hcfv8&dl=0

In this link you can see the video that describes the problem that i see.

[Torsten]

And you don't think that a steady-state could be reached, i.e. an equilibrium between growth and death of the algae or whatever process is modelled by your equations ?

[Alfonso]

Yes, it could. But im changing also the initial parameters P, N and C and always it stops around z=4

[Alfonso]

But what i wanna try i also to implement the equations (9) in the file pdf algae and try to see what happens. Can you help me to do it?

[Torsten]

f1 = some function of P(j)
f2 = some function of N(j)
f3 = some function of C(j)
N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j) + S_N)
P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j) + S_P)
C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j) + S_C)

I think within the formula

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));

integral_term has has to be replaced by the integrand f1*f2*f3*g*w itself and does not enter as integral(f1*f2*f3*g*w dz) ? Further, why is "integral_term" multiplied by w(2:nz,j) ?

Should simply be

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)

I guess.

[Alfonso]

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 = 8000; % [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;
%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

    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) = 0%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(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
%plot(z,flipud(w(:,end)))
end

Done! But at T=1000 i have problems, after this time you see that the code is instable

[Torsten]

P = 3;
N = 4;
C = 3 ;

Are these really the correct values for P, N and C at time t = 0 or are these different parameters ?

If they are correct, please use them as

P = zeros(nt+1,1)
N = zeros(nt+1,1)
C = zeros(nt+1,1)
P(1) = 3;
N(1) = 4;
C(1) = 3;

before you enter the for-loop.

And as suggested by me and other forum members several times now: choose a smaller time step (thus enlarge nt). Then the stability problems will vanish.

And to avoid the flipping of the solution for w, simply set

w = zeros(nz+1,nt+1);
w(1,:) = 5;

when initializing w and remove the lines

    w(1,j+1) = 0%w(2,j+1);
    w(nz+1,j+1) = 5;%w(nz,j+1); %Concentration at z=0   

within the for-loop.

The condition w(z=Inf) = 0 is not natural - you should take care that the value for w at z = Z never influences what happens within 0 <= z < Z. A condition dw/dz = 0 at z = Z would make more sense.