Im working with algae growth linked to the light. But i have some problems on my code
1 次查看(过去 30 天)
显示 更早的评论
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1545522/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1545527/image.png)
Do you have some suggestions to solve my ploblems? Thank you very much
As attachment you can find the file that im following.
1 个评论
回答(1 个)
Ganesh
2023-12-29
The issue you are facing is due to the following lines:
% At T=1001, the loop conditions are met and results into an infinite loop
while (2*r) > st
nt = nt+1;
end
The while statement has an update condition that doesn't seem to modify the loop condition. I understand you may need to calculate new values of `r` based on the updated value of nt, for which you may need to use a make-do "do-while loop".
Kindly refer to the Answer to know more on implementing the same: https://in.mathworks.com/matlabcentral/answers/115403-do-while-loop-in-matlab
While this serves as the answer to your question, it would be a better option for you to learn the debugging techniques as mentioned by @Image Analyst.
Kindly refer to the following documentation to know more on debugging MATLAB Code Files:
Hope this helps!
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!