how to use derivative of function using gradient?

31 views (last 30 days)
hi, im a student trying to solve a mathematical problem using MATLAB, the code consists ODE's (ordinary differential equations). the code is using a numerical analysis in order to solve all the ODE's simoultenously. One of the ODE's using a derivative of a parameter called "par". the parameter's value is changing along z axis (the problem defined only for z axis). this is the line:
dydz(1) = -(gradient(par,z))*(y(1)/par);
unfortenately when i display all the values of "par" using: disp(par); i get a lot of different numbers as expected. but when i use disp(gradient(par,z)); i get only zero's 0.
what am i doing wrong? is gradient function gets the gradient of a single point each time and returns zero? how do i fix this?
Thanks !
  1 Comment
Gal Shaked
Gal Shaked on 21 Sep 2022
i attached the code, i erased many lines so you can focus on the relevant lines.
the problem is in line: dydz(1) = -(gradient(dens_mix,z))*(y(1)/dens_mix);
while the parameter is from: dens_mix= partial_dens(1)+partial_dens(2)+partial_dens(3)+partial_dens(4);
and:
partial_dens = [0 0 0 0 ];
for i = 1:length(partial_dens)
partial_dens(i) = (partial_press(i)*molecular_weights(i))/(Rconstant*y(11));
end
the parameters y(11) and partial_press(i) -values are changing along z axis while the others are constants

Sign in to comment.

Answers (2)

Torsten
Torsten on 21 Sep 2022
Edited: Torsten on 21 Sep 2022
Add the variable "dens_mix" as an algebraic equation to your system (e.g. as y(13)):
dydz(13) = y(13) - (partial_dens(1)+partial_dens(2)+partial_dens(3)+partial_dens(4));
and supply the mass matrix of your ODE system as
function M = Mass(t,y)
M = eye(13);
M(13,13) = 0;
M(1,13) = y(1)/y(13);
end
If follows that in the vector dydz that you prescribe in ODEBVP, you have to set
dydz(1) = 0
If you don't use an ODE integrator (like ODE45 or ODE15S), but a BVP solver (like BVP4C), come back here. In this case, a different solution will be necessary.
  15 Comments
Gal Shaked
Gal Shaked on 26 Sep 2022
Edited: Gal Shaked on 26 Sep 2022
no, i wrote dydz(13). it came frmo overall mass balance:
(rho means density- dens_mix), and u is y(1) and rho is y(13)
after differentiation we get: - that equation definens dydz(1) equation and dydz(13) equation.
how do i use that toolbox in order to define dens_mix as a function? it's value is changing each time the code is actually trying to solve the ODE system, and it is just a sum of 4 parameters:
dens_mix= partial_dens(1)+partial_dens(2)+partial_dens(3)+partial_dens(4)
and thank again, really appreciate that.

Sign in to comment.


Torsten
Torsten on 26 Sep 2022
Edited: Torsten on 26 Sep 2022
Forget about y(13) and use
dydz(2) = y(3); %Specific mass balance (CH4)
dydz(3) = ((y(1)*y(3))+(dens_cat_weighted*r_CH4))/(D_m+D_l); %Specific mass balance (CH4)
dydz(4) = y(5); %Specific mass balance(H2O)
dydz(5) = ((y(1)*y(5))+(dens_cat_weighted*r_H2O))/(D_m+D_l); %Specific mass balance(H2O)
dydz(6) = y(7); %Specific mass balance(CO2)
dydz(7) = ((y(1)*y(7))+(dens_cat_weighted*r_CO2))/(D_m+D_l); %Specific mass balance(CO2)
dydz(8) = y(9); %Specific mass balance (H2)
dydz(9) = ((y(1)*y(9))+(dens_cat_weighted*r_H2))/(D_m+D_l); %Specific mass balance(H2)
dydz(10) = -(1/1000)*(((dyn_visc_mix/K_DPC)*y(1))+((dens_mix/K_nDC)*(y(1)^2))); % Momentum balance
dydz(11) = y(12); %energy balance (temprature)
dydz(12) = ((enthalpy*r*dens_cat_weighted*1000) - (flux/H) + (specific_heat*dens_mix*y(1)*y(12)))/(k_m); %energy balance
d_dens_mix_dz = (dydz(10)*y(11)-y(10)*dydz(11))/(Rconstant*y(11)^2)*(...
molecular_weights(1)*y(2)/(y(2)+y(4)+y(6)+y(8))+...
molecular_weights(2)*y(4)/(y(2)+y(4)+y(6)+y(8))+...
molecular_weights(3)*y(6)/(y(2)+y(4)+y(6)+y(8))+...
molecular_weights(4)*y(8)/(y(2)+y(4)+y(6)+y(8)))+...
(y(10)/(Rconstant*y(11))*(...
molecular_weights(1)*((dydz(2)*(y(2)+y(4)+y(6)+y(8))-...
y(2)*(dydz(2)+dydz(4)+dydz(6)+dydz(8)))/...
(y(2)+y(4)+y(6)+y(8))^2)+...
molecular_weights(2)*((dydz(4)*(y(2)+y(4)+y(6)+y(8))-...
y(4)*(dydz(2)+dydz(4)+dydz(6)+dydz(8)))/...
(y(2)+y(4)+y(6)+y(8))^2)+...
molecular_weights(3)*((dydz(6)*(y(2)+y(4)+y(6)+y(8))-...
y(6)*(dydz(2)+dydz(4)+dydz(6)+dydz(8)))/...
(y(2)+y(4)+y(6)+y(8))^2)+...
molecular_weights(4)*((dydz(8)*(y(2)+y(4)+y(6)+y(8))-...
y(8)*(dydz(2)+dydz(4)+dydz(6)+dydz(8)))/...
(y(2)+y(4)+y(6)+y(8))^2));
dydz(1) = -d_dens_mix_dz*(y(1)/dens_mix); %Overall mass balance
You might want to check the derivative using symbolic computation:
syms z y10(z) y11(z) y2(z) y4(z) y6(z) y8(z) M1 M2 M3 M4 Rconstant
dens_mix = y10/(Rconstant*y11)*(M1*y2+M2*y4+M3*y6+M4*y8)/(y2+y4+y6+y8);
d_dens_mix_dz = simplify(diff(dens_mix,z))
d_dens_mix_dz(z) = 
  6 Comments
Gal Shaked
Gal Shaked on 28 Sep 2022
  1. ok ill try this
  2. ok got it
  3. yes sorry, i mean i used what you wrote and the code didnt respond for a while until i got an error: "Unable to solve the collocation equations -- a singular Jacobian encountered."
thanks again

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!