i am getting the following issue when i am trying to calculate gradient and divergence in matlab' ode45.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%this is 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.1; % Plz confirm this
IC=[q_d1_0 q_d2_0 q_r_0 c_small_0]; % IC is the set of all the Initial conditions
%% Function call
%*****************************
theta_values = [80+273, 100+273, 120+273]; % Define different theta values
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
%% Plot cdot vs. time for different temperatures
% Initialize a cell array to store cdot and k_small data for each temperature
cdot_data = cell(1, length(theta_values));
k_small_data = cell(1, length(theta_values)); 
for i = 1:length(theta_values)
    theta = theta_values(i);
    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
    % Solve the ODE for the current theta
    [IVsol, DVsol] = ode45(@(IV, DV) FUNC_DEF(IV, DV, theta), t_span, IC);
    
    % Calculate k_small and der_c_small_def for the current theta
    k_small = (k01*(1-((DVsol(:, 1)+DVsol(:, 2))/2))+k02*(1-DVsol(:, 3)))*exp(-(E_k / (R * theta)));
    der_c_small_def = ((-k_small.*DVsol(:, 4))/rho_0);
    
    % % Plot cdot vs. time for the current temperature
    % semilogx(IVsol, der_c_small_def, 'LineWidth', 2, 'DisplayName', [num2str(theta_values(i) - 273), '°C']);
    
    % Store cdot and k_small data for this theta
    cdot_data{i} = der_c_small_def;
    k_small_data{i} = k_small;
    % Create a new figure for each temperature
    % ********************************
    figure();  
    semilogx(IVsol, der_c_small_def, 'LineWidth', 2);
    xlabel("t [s]", 'FontSize', 12,'FontWeight','bold');
    ylabel("$\dot{c}$ [mol/kg/s]", 'Interpreter', 'latex', 'FontSize', 12,'FontWeight','bold');
    set(gca, 'XScale', 'log');
    grid on;
    title(['Rate of Oxygen Absorption vs. Time at ', num2str(theta - 273), '°C'], 'FontSize', 14);
    ax = gca;%Set axis font size
    ax.XAxis.FontSize = 24;
    ax.YAxis.FontSize = 24;
    figure();  
    semilogx(IVsol, DVsol(:,4), 'LineWidth', 2);
    xlabel("t [s]", 'FontSize', 12,'FontWeight','bold');
    ylabel("c", 'Interpreter', 'latex', 'FontSize', 12,'FontWeight','bold');
    set(gca, 'XScale', 'log');
    grid on;
    title(['Oxygen concentration vs. Time at ', num2str(theta - 273), '°C'], 'FontSize', 14);
    ax = gca;%Set axis font size
    ax.XAxis.FontSize = 24;
    ax.YAxis.FontSize = 24;
    
    figure();  
    semilogx(IVsol, k_small, 'LineWidth', 2);
    xlabel("t [s]", 'FontSize', 12,'FontWeight','bold');
    ylabel("k", 'Interpreter', 'latex', 'FontSize', 12,'FontWeight','bold');
    set(gca, 'XScale', 'log');
    grid on;
    title(['Reaction rate vs. Time at ', num2str(theta - 273), '°C'], 'FontSize', 14);
    ax = gca;%Set axis font size
    ax.XAxis.FontSize = 24;
    ax.YAxis.FontSize = 24;
    % ********************************
    % Save the results for this temperature into a MAT file
    results = struct('IVsol', IVsol, 'DVsol', DVsol, 'cdot', der_c_small_def, 'k_small', k_small);
    filename = ['Diffusion_specific_results_', num2str(theta - 273), 'C.mat'];
    save(filename, 'results');
end

Error using divergence
Not enough input arguments.

Error in solution>FUNC_DEF (line 177)
div_dc_grad_c = divergence(d_c * grad_c_small);

Error in solution>@(IV,DV)FUNC_DEF(IV,DV,theta) (line 44)
[IVsol, DVsol] = ode45(@(IV, DV) FUNC_DEF(IV, DV, theta), t_span, IC);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% Create one figure for all temperatures : plz use lines 53-54 mandatorily as well
% ********************************
% xlabel("t [s]", 'FontSize', 12,'FontWeight','bold');
% ylabel("$\dot{c}$ [mol/kg/s]", 'Interpreter', 'latex', 'FontSize', 12,'FontWeight','bold');
% set(gca, 'XScale', 'log');
% grid on;
% legend('Location', 'northeast');
% title("Rate of Oxygen Absorption vs. Time for Different Temperatures", 'FontSize', 14);
% ax = gca;%Set axis font size
% ax.XAxis.FontSize = 24;
% ax.YAxis.FontSize = 24;
% ********************************  
%% Stop timing and display elapsed time
%*****************************
elapsed_time = toc;
fprintf('Total time taken: %.4f seconds\n', elapsed_time);
%plot for q_d1, q_d2 ,q_r ,c_small  
figure()
plot(IVsol, DVsol), 
grid on
xlabel('IV'), ylabel('DV')
legend({'q_d1','q_d2','q_r','c_small'}, 'location', 'east')
%plot for q_d, q_r ,c_small ,k_small  
q_d_new = DVsol(:,1)+DVsol(:,2);
figure()
plot(IVsol, q_d_new,'*g' ,IVsol, DVsol(:,3),'-b' ,IVsol, DVsol(:,4),'--r') 
legend({'q_d','q_r','c_small'}, 'location', 'east')
%curve fitting plot for c
DVsol_cplot=DVsol(:,4);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%this is FUNC_DEF : (the function defination)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dDVdIV] = FUNC_DEF(IV, DV, theta)
% I : Independent Var, D : Dependent Var; dDVdt_DEF: funnction defination
% IV , I , IVsol : Independent Variables
% DV , D , DVsol : Independent Variables
% Defining constants
%*****************************
% theta = 80+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_d)+k02*(1-q_r))*exp(-(E_k / (R * theta))); % expr of k_small
grad_c_small = gradient(c_small);
div_dc_grad_c = divergence(d_c * grad_c_small);
der_c_small_def = ((div_dc_grad_c - k_small .* c_small) / rho_0);% DE of c_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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

i am getting this error:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Error using divergence

Not enough input arguments.

Error in FUNC_DEF (line 49)

div_dc_grad_c = divergence(d_c * grad_c_small);

Error in FUNC_MAIN>@(IV,DV)FUNC_DEF(IV,DV,theta) (line 48)

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

Error in odearguments (line 92)

f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in FUNC_MAIN (line 48)

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

[Star Strider]

the divergence function requires least 2 matrices for its arguments that are at least each a (2x2) matrix. You are giving it one scalar (equal to 0).

Answer 1

[Ganesh]

Hello @S Ashish,

The divergence() function calculates the divergance by taking two matrices for a 2D space and 3 matrices for a 3D space. At the very least, 2 arguments need to be passed to the divergence function. In the line

divergence(d_c * grad_c_small);

you pass only one argument, which gives rise to the issue.

Please do also note that the grad_c_small is 0, as the gradient of a constant value(c_small) is 0. This further makes the passed argument a scalar and not a matrix. Kindly check the implementation of your code to resolve this.

Please find linked, the documentation of gradient and divergance functions.

Thanks and Regards,

Ganesh Saravanan