i am getting the following issue when i am trying to calculate gradient and divergence in matlab' ode45.
2 views (last 30 days)
Show older comments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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
% 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);
1 Comment
Star Strider
on 5 Sep 2023
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).
Answers (1)
Ganesh
on 12 Sep 2023
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.
Thanks and Regards,
Ganesh Saravanan
0 Comments
See Also
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!