MATLAB: Dsolve – Warning: Unable to find symbolic solution.

dsolveMATLABsymbolic solution

% 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?

Best Answer

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
Related Question