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

Question

S Ashish 2023-8-25

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2012722-i-am-using-ode45-to-solve-a-system-of-4-odes-my-dvsol-that-contains-all-results-for-my-dpendent-v

评论： Sam Chak 2023-8-29

采纳的回答： Sam Chak

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Torsten 2023-8-25

编辑：Torsten 2023-8-25

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 2023-8-29

编辑：S Ashish 2023-8-29

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Sam Chak 2023-8-25

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2012722-i-am-using-ode45-to-solve-a-system-of-4-odes-my-dvsol-that-contains-all-results-for-my-dpendent-v#answer_1293877

编辑：Sam Chak 2023-8-25

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

S Ashish 2023-8-29

移动：Sam Chak 2023-8-29

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

Sam Chak 2023-8-29

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!

请先登录，再进行评论。