MATLAB Answers

dsolve - Warning: Unable to find symbolic solution.

28 views (last 30 days)
% Parameters
mu_d = 0; % specific decay rate
mu_N = 2.69; % maximum specific nitrate uptake rate
K_N = 0.8; % half velocity coefficient
k_q = 19.6; % minimum nitrogen quota
theta = 6.69; % modified kinetic constant for biolipid syntheiss
epsilon = 0.001; % modified parameter taking into account the complex effect of lipid esterification
gamma = 7.53; % modified kinetic constant for biolipid consumption
tau = 0.138; % dimensionless kinetic parameter
delta = 9.9; % dimensionless kinetic parameter
phi = -0.456; % dimensionless kinetic parameter
I_0 = 80; % incident light intensity micromol/m^2.s from experimental data
alpha = 196.4; % cell absorption coefficient
beta = 0; % bubble scattering coefficient
L = 0.044; % reactor length
k_i = 100; % light inhibition terms
mu_max = 0.36; % maximum specific growth rate in 1 h
k_s = 91.2; % light saturation term in micromol/m^2.s
syms X(t) N(t) q(t) f(t)
% Modified Beer-Lambert law
I = sym(zeros(11,1));
for z = 0:10
I(z+1) = I_0 * exp(-(alpha*X + beta) * (z*L)/10);
end
% Trapezoidal rule
mu_m_first = I(1)/(I(1)+k_s+I(1)^2/k_i);
mu_m_last = I(11)/(I(11)+k_s+I(11)^2/k_i);
mu_m_middle = 0;
for n = 2:10
mu_m_middle = mu_m_middle + (I(2)/(I(2)+k_s+I(2)^2/k_i));
end
mu_m = mu_max/20 * (mu_m_first + 2*mu_m_middle + mu_m_last);
% Specific growth rate
mu_0 = mu_m * (1 - k_q/q);
% Algal biomass growth rate
ode1 = diff(X) == mu_0*X - mu_d*X;
% Nitrate consumption
ode2 = diff(N) == -mu_N * (N/(N+K_N)) * X;
% Accumulation rate of nitrogen quota
ode3 = diff(q) == mu_N * (N/(N+K_N)) - mu_m * (1-k_q/q) * q;
% FAME production
ode4 = diff(f) == mu_m * (theta*q - epsilon*f) * (1-k_q/q) - gamma*mu_N*(N/(N+K_N));
odes = [ode1; ode2; ode3; ode4];
% Chlorophyll fluorescence
Y_II = exp(tau*q)/(exp(tau*q)+delta) + phi;
% Initial conditions
cond1 = X(0) == 0.178333333; % initial biomass concentration in g/L from experimental data
cond2 = N(0) == 35.01666667; % initial culture nitrate concentration in mg/L from experimental data
cond3 = q(0) == 79.85; % initial nitrogen quota in mg/g
cond4 = f(0) == 120.3100949; % initial FAME yield in %wt from experimental data
conds = [cond1; cond2; cond3; cond4];
I get this warning if I use this line.
% Warning: Unable to find symbolic solution.
S = dsolve(odes, conds);
If I use the following line instead, I get an extra error.
% Warning: Unable to find symbolic solution.
% Error using sym/subsindex (line 853)
% Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be symbolic variables, and function body must be sym expression.
[XSol(t), NSol(t), qSol(t), fSol(t)] = dsolve(odes, conds);
% plot
fplot(XSol)
hold on
fplot(NSol)
fplot(qSol)
fplot(fSol)
grid on
legend('XSol', 'NSol', 'qSol', 'fSol', 'Location','best')
What have I done wrong here?

  0 Comments

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 20 Nov 2020
Why not just solve them numerically:
% Parameters
tau = 0.138; % dimensionless kinetic parameter
delta = 9.9; % dimensionless kinetic parameter
phi = -0.456; % dimensionless kinetic parameter
% Initial conditions
X(1) = 0.178333333; % initial biomass concentration in g/L from experimental data
N(1) = 35.01666667; % initial culture nitrate concentration in mg/L from experimental data
q(1) = 79.85; % initial nitrogen quota in mg/g
f(1) = 120.3100949; % initial FAME yield in %wt from experimental data
B0 = [X(1) N(1) q(1) f(1)];
tspan = [0 200];
[t, B] = ode45(@odefn, tspan, B0);
X = B(:,1); N = B(:,2); q = B(:,3); f = B(:,4);
% Chlorophyll fluorescence
Y_II = exp(tau*q)./(exp(tau*q)+delta) + phi;
subplot(2,1,1)
plot(t,X),grid
xlabel('t'),ylabel('Algal biomass growth rate')
subplot(2,1,2)
plot(t,N),grid
xlabel('t'),ylabel('Nitrate consumption')
figure
subplot(2,1,1)
plot(t,q),grid
xlabel('t'),ylabel('Accum. rate of nitrogen quota')
subplot(2,1,2)
plot(t,f),grid
xlabel('t'),ylabel('FAME production')
figure
plot(t,Y_II),grid
xlabel('t'),ylabel('Chlorophyll fluorescence')
function dBdt = odefn(~,B)
% Data
mu_d = 0; % specific decay rate
mu_N = 2.69; % maximum specific nitrate uptake rate
K_N = 0.8; % half velocity coefficient
k_q = 19.6; % minimum nitrogen quota
theta = 6.69; % modified kinetic constant for biolipid syntheiss
epsilon = 0.001; % modified parameter taking into account the complex effect of lipid esterification
gamma = 7.53; % modified kinetic constant for biolipid consumption
I_0 = 80; % incident light intensity micromol/m^2.s from experimental data
alpha = 196.4; % cell absorption coefficient
beta = 0; % bubble scattering coefficient
L = 0.044; % reactor length
k_i = 100; % light inhibition terms
mu_max = 0.36; % maximum specific growth rate in 1 h
k_s = 91.2; % light saturation term in micromol/m^2.s
% Extract parameters
X = B(1); N = B(2); q = B(3); f = B(4);
% Modified Beer-Lambert law
I = zeros(11,1);
for z = 0:10
I(z+1) = I_0 * exp(-(alpha*X + beta) * (z*L)/10);
end
% Trapezoidal rule
mu_m_first = I(1)/(I(1)+k_s+I(1)^2/k_i);
mu_m_last = I(11)/(I(11)+k_s+I(11)^2/k_i);
mu_m_middle = 0;
for n = 2:10
mu_m_middle = mu_m_middle + (I(2)/(I(2)+k_s+I(2)^2/k_i));
end
mu_m = mu_max/20 * (mu_m_first + 2*mu_m_middle + mu_m_last);
% Specific growth rate
mu_0 = mu_m * (1 - k_q/q);
% Rate equations
dBdt = [mu_0*X - mu_d*X;
-mu_N * (N/(N+K_N)) * X;
mu_N * (N/(N+K_N)) - mu_m * (1-k_q/q) * q;
mu_m * (theta*q - epsilon*f) * (1-k_q/q) - gamma*mu_N*(N/(N+K_N))];
end

  3 Comments

StillANovice
StillANovice on 20 Nov 2020
Hi Alan, thanks for this!
But I am still getting the same error:
% Error using sym/subsasgn (line 959)
% Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be symbolic
% variables, and function body must be sym expression.
% Error in chbe_6500 (line 7)
% X(1) = 0.178333333; % initial biomass concentration in g/L from experimental data
Alan Stevens
Alan Stevens on 20 Nov 2020
The code I posted doesn't use symbolic commands at all! Just copy and paste into a completely blank new script. Here's one of the graphs that results:
StillANovice
StillANovice on 20 Nov 2020
sorry, I forgot to clear up my workspace!
Thanks Alan!

Sign in to comment.

More Answers (0)

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!