i am getting the following issue when i am trying to calculate gradient and divergence in matlab' ode45.
6 次查看(过去 30 天)
显示 更早的评论
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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 个评论
Star Strider
2023-9-5
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).
回答(1 个)
Ganesh
2023-9-12
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Annotations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!