I am using ODE45 to solve a system of 4 ODEs. my DVsol ( that contains all results for my Dpendent Variables) are showing zero.

This is the FUNC_MAIN ( the Function that calls FUNC_DEF) :

clc
clear all  
close all 

% Start timing
tic
%% Inputs
%*****************************
t_span=[0 1e7]; %time span
q_d1_0=0; % Initial conditions for q_d1,q_d2, q_r and c_small
q_d2_0=0;
q_r_0=0;
c_small_0=0;
IC=[q_d1_0 q_d2_0 q_r_0 c_small_0]; % IC is the set of all the Initial conditions

%% Function call
%*****************************
% %auto time step

% [IVsol, DVsol] = ode45(@(IV, DV) FUNC_DEF(IV, DV), t_span, IC);
% [IVsol, DVsol] = ode45('FUNC_DEF', t_span, IC);


% adjusted time step
%options = odeset('MaxStep', 1); % Adjust MaxStep as needed
[IVsol, DVsol] = ode15s(@(IV, DV) FUNC_DEF(IV, DV), t_span, IC);%, options);  

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

1.0e-19 * 0.9513 0.1362 0.0024 -0.0339

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

dDVdIV = 4×1

0 0 0 0

plot(IVsol,DVsol)

progress = (IVsol(end) - t_span(1)) / (t_span(2) - t_span(1)) * 100; % Display progress
fprintf('Integration progress: %.2f%%\n', progress);

Integration progress: 100.00%

%% Save results in a MAT file
%*****************************
results = struct('IVsol', IVsol, 'DVsol', DVsol);
%save('Diffusion_specific_results_120C.mat', 'results');
% save("Diffusion_specific_results.mat")  
% disp('Results saved');

%% Stop timing and display elapsed time
%*****************************
elapsed_time = toc;
fprintf('Total time taken: %.4f seconds\n', elapsed_time);

Total time taken: 0.5132 seconds

This is the FUNC_DEF ( the Function defination program) :

function [dDVdIV] = FUNC_DEF(IV, DV)  
% i : Independent Var, D : Dependent Var; dDVdt_DEF: funnction defination

% IV , I , IVsol : Independent Variables
% DV , D , DVsol : Independent Variables
% Defining constants
%*****************************
theta = 120+273;                        % absloute temperature in Kelvin
k01=2.776e7;                            % in /sec
k02=2.136e10;                           % in /sec
E_k=1.0513e5;                           % in J/mol
% c_0 = 1;                                % dimensionless local oxygen concentration: ASSUMED
R = 8.314;                              % universal gas constant in J/mol
rho_0=1000;                             % Density of NBR in the reference configuration in kg/m^3
d_c_0=1.1e-10;                          % in m^2/sec
gamma_c=72.196;
nu_d1 = 5.3955e6;                        % in /sec
E_d1 = 8.9732e4;                         % in J/mol
nu_d2 = 5.401e7;                        % in /sec
E_d2 = 1.0361e5; 
cons_d1_wo_c = nu_d1 * exp(-(E_d1 / (R * theta)));   % cons_d1_wo_c=cons_d1 without c
cons_d2_wo_c = nu_d2 * exp(-(E_d2 / (R * theta)));   % cons_d2_wo_c=cons_d2 without c
nu_r = 8.745e4;                          % in /sec 
E_r = 9.5852e4;                          % in J/mol
cons_r_wo_c = nu_r * exp(-(E_r / (R * theta)));      % cons_r_wo_c=cons_r without c
%*****************************
% Defining Dependent variables
%*****************************
q_d1=DV(1);
q_d2=DV(2);
q_r=DV(3);
c_small=DV(4); 
%*****************************

% Calculate the constants including c
%*****************************
cons_d1_with_c = c_small .* cons_d1_wo_c ;
cons_d2_with_c = c_small.* cons_d2_wo_c ;
cons_r_with_c = c_small.* cons_r_wo_c ;
%*****************************
% Calculate the derivatives of DVs and Extra variables defination
%*****************************
der_q_d1_def=cons_d1_with_c .* (1 - q_d1); % DE of q_d1
der_q_d2_def=cons_d2_with_c .* (1 - q_d2); % DE of q_d2
q_d=(q_d1+q_d2)/2;                         % expr of q_d
der_q_r_def=cons_r_with_c .* (1 - q_r);    % DE of q_r
d_c = d_c_0*exp(-(gamma_c*q_r));           % expr of d_c
k_small=(k01*(1-((q_d1+q_d2)/2))+k02*(1-q_r))*exp(-(E_k / (R * theta))); % expr of k_small
der_c_small_def=((-k_small.*c_small)/rho_0); % expression for DE of c_dot_small
%***************************** 
% Defining Function defination
%*****************************
dDVdIV=[der_q_d1_def;der_q_d2_def;der_q_r_def;der_c_small_def]
% dDVdIV=@(DV,IV)[der_q_d1_def;der_q_d2_def;der_q_r_def;der_c_small_def];
%*****************************
end

[Torsten]

As you can see, the derivatives of your solution variables all turn out to be 0. So starting with 0 and adding 0 gives 0.

[S Ashish]

Yes, you are right. Setting a non-zero IC did the trick . Thank You.

Answer 1

[Sam Chak]

Hi @S Ashish

Update: In short, the zeros are the equilibrium points of the system. Here are Trajectories obtained from the non-zero initial values.

%% Inputs
%*****************************
t_span    = [0 1e7];    % time span
q_d1_0    =  50;        % Initial conditions for q_d1,q_d2, q_r and c_small
q_d2_0    = -20;
q_r_0     =  10;
c_small_0 =   4;
IC        = [q_d1_0 q_d2_0 q_r_0 c_small_0]; % IC is the set of all the Initial conditions

%% Function call
%*****************************
% %auto time step

% [IVsol, DVsol] = ode45(@(IV, DV) FUNC_DEF(IV, DV), t_span, IC);
% [IVsol, DVsol] = ode45('FUNC_DEF', t_span, IC);


% adjusted time step
%options = odeset('MaxStep', 1); % Adjust MaxStep as needed
[IVsol, DVsol] = ode15s(@(IV, DV) FUNC_DEF(IV, DV), t_span, IC);%, options);  

plot(IVsol, DVsol), grid on, legend('show', 'location', 'east')
xlabel('IV'), ylabel('DV')

function [dDVdIV] = FUNC_DEF(IV, DV)  
    % i : Independent Var, D : Dependent Var; dDVdt_DEF: funnction defination
    
    % IV , I , IVsol : Independent Variables
    % DV , D , DVsol : Independent Variables
    % Defining constants
    %*****************************
    theta   = 120+273;      % absloute temperature in Kelvin
    k01     = 2.776e7;      % in /sec
    k02     = 2.136e10;     % in /sec
    E_k     = 1.0513e5;     % in J/mol
    % c_0   = 1;            % dimensionless local oxygen concentration: ASSUMED
    R       = 8.314;        % universal gas constant in J/mol
    rho_0   = 1000;         % Density of NBR in the reference configuration in kg/m^3
    d_c_0   = 1.1e-10;      % in m^2/sec
    gamma_c = 72.196;
    nu_d1   = 5.3955e6;     % in /sec
    E_d1    = 8.9732e4;     % in J/mol
    nu_d2   = 5.401e7;      % in /sec
    E_d2    = 1.0361e5; 
    nu_r    = 8.745e4;      % in /sec 
    E_r     = 9.5852e4;     % in J/mol
    
    cons_d1_wo_c = nu_d1 * exp(-(E_d1 / (R * theta)));   % cons_d1_wo_c=cons_d1 without c
    cons_d2_wo_c = nu_d2 * exp(-(E_d2 / (R * theta)));   % cons_d2_wo_c=cons_d2 without c
    cons_r_wo_c  = nu_r  * exp(-(E_r  / (R * theta)));   % cons_r_wo_c=cons_r without c
    
    %*****************************
    % Defining Dependent variables
    %*****************************
    q_d1    = DV(1);
    q_d2    = DV(2);
    q_r     = DV(3);
    c_small = DV(4); 
    
    %*****************************
    % Calculate the constants including c
    %*****************************
    cons_d1_with_c = c_small.*cons_d1_wo_c;
    cons_d2_with_c = c_small.*cons_d2_wo_c;
    cons_r_with_c  = c_small.*cons_r_wo_c;
    
    %*****************************
    % Calculate the derivatives of DVs and Extra variables defination
    %*****************************
    der_q_d1_def    = cons_d1_with_c .* (1 - q_d1); % DE of q_d1
    der_q_d2_def    = cons_d2_with_c .* (1 - q_d2); % DE of q_d2
    q_d             = (q_d1 + q_d2)/2;                         % expr of q_d
    der_q_r_def     = cons_r_with_c .* (1 - q_r);    % DE of q_r
    d_c             = d_c_0*exp(-(gamma_c*q_r));           % expr of d_c
    k_small         = (k01*(1-((q_d1+q_d2)/2))+k02*(1-q_r))*exp(-(E_k / (R * theta))); % expr of k_small
    der_c_small_def = ((-k_small.*c_small)/rho_0); % expression for DE of c_dot_small
    
    %***************************** 
    % Defining Function defination
    %*****************************
    dDVdIV = [der_q_d1_def;
              der_q_d2_def;
              der_q_r_def;
              der_c_small_def];
    
    % dDVdIV=@(DV,IV)[der_q_d1_def;der_q_d2_def;der_q_r_def;der_c_small_def];
    %*****************************
end

[S Ashish]

Yes, you are right. Setting a non-zero IC did the trick . Thank You.

[Sam Chak]

You are welcome, @S Ashish. If you find the solution helpful, please consider clicking 'Accept' ✔ the answer to close the issue when the problem is resolved. Thanks a bunch!