how to use derivative of function using gradient?

29 次查看(过去 30 天)
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 个评论
Gal Shaked
Gal Shaked 2022-9-21
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

请先登录,再进行评论。

回答(2 个)

Torsten
Torsten 2022-9-21
编辑:Torsten 2022-9-21
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 个评论
Torsten
Torsten 2022-9-25
编辑:Torsten 2022-9-25
The use of the symbolic toolbox is only meant to differentiate the expression for dens_mix with respect to z. Once you have this derivative, you can copy it in your existing code.
But you write above
dydz(13) = something
So you already differentiated y(13) = dens_mix with respect to z and the result is the right-hand side ?
I think you wrote y(13) and not dydz(13), didn't you ?
Gal Shaked
Gal Shaked 2022-9-26
编辑:Gal Shaked 2022-9-26
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.

请先登录,再进行评论。


Torsten
Torsten 2022-9-26
编辑:Torsten 2022-9-26
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 个评论
Torsten
Torsten 2022-9-28
编辑:Torsten 2022-9-28
how can i use a variable from one MATLAB file in another one? as in the picture, a variable from ODEBVP.m that i want to use in bvp4bc.m
I don't know how the two functions are connected. If nothing helps, use a global variable.
after i solved the velocity u(z), how can i present it in a graph? how do i returning the u(z) as a funciton result? i am presenting the other solved parameters (calculated in ODEBVP.m) in the main file BVP.m just as written below but im not sure how to return the u(z)
You must recalculate density_mix(z) for all z-positions from the other solution variables sol.y and then calculate u from u(z) = (u*density_mix)(@z=0)/density_mix(z)
hi, i know you gave me the 13th equation, the thing is- the code can't run it, it getting stuck by errors.
I didn't give you the 13th equation. I just gave you the formula how to replace -(gradient(dens_mix,z)) in the definition of dydz(1). A 13th equation is not necessary.
Gal Shaked
Gal Shaked 2022-9-28
  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

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by