how to use derivative of function using gradient?

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 !

[Gal Shaked]

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

Answer 1

[Torsten]

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.

[Gal Shaked]

Thank you for your comment, im sorry i didnt mention this before but im using BVP solver- i think it is actually BVP4C

im not sure how to verify which kind of BVP solver im specifically using, but i think it is BVP4C.

[Torsten]

First define dy(2)/dz,...,dy(12)/dz as you did in your code.

Then you will have to differentiate using pencil and paper the identity

dens_mix = (partial_dens(1)+partial_dens(2)+partial_dens(3)+partial_dens(4));

with respect to z (and thus express d(dens_mix)/dz in terms of the z-derivatives of the variables from y(1),...,y(12) that are involved).

If d(dens_mix)/dz comes out to be some function f(dy(1)/dz,dy(2)/dz,...,dy(12)/dz), you can insert this expression in the first ODE:

dy(1)/dz = -f(dy(1)/dz,dy(2)/dz,...,dy(12)/dz)*y1/dens_mix

Since possibly both the left and the right-hand side contain dy1/dz, it might be necessary to solve for it.

[Gal Shaked]

thank you for your help, but i dont really understand how to define: d(dens_mix)/dz, its equation is including the same parameters as dydz(1). i tried something but got a constant solution to dens_mix which means something wrong

[Torsten]

What comes out if you express d(dens_mix)/dz using y(1),...,y(12),dydz(1),...,dydz(12) ?

[Gal Shaked]

the code solved 'dens_mix' as a constant value, equals to the value i defind as a initial boundary condition.

i assume that it relates to the fact that both equations 1 and 13 (13 is the equation represents dens_mix) are mathematically the same, just rearranged.

[Torsten]

What is equation 13 ? I do not yet know your new code.

[Gal Shaked]

i attached a new simplyfied version of a part of the code, it's not for running but only to get the idea (i erased many lines that are irrelevant for the problem)

[Torsten]

We need dens_mix written in y(1),...,y(12). Can you establish this ?

I would have deduced it from your code, but it's missing how molar_fractions are derived from the solution variables.

[Gal Shaked]

i believe i can rephrase dydz(13):

dydz(13) = (((dydz(1))/(Rconstant*y(11)*y(1)))*(molecular_weights(1)*partial_press(1)+ ...
    molecular_weights(2)*partial_press(2)+molecular_weights(3)*partial_press(3)+ ...
    molecular_weights(4)*partial_press(4)));

• partial_pressure - value that changes according to z
• y(11) - temperature changes according to z
• y(1) velocity changes according to z
• Rconstant and molecular_weights are constants

but still i get a constant velocity

would it help if i will post or send you privetly the full code?

thanks again

[Torsten]

Express density_mix in terms of y(1),..,y(12) and constants/parameters.

partial_press depends on y(10) and molar_fractions where (I guess) molar_fractions also depends on y(1),...y(12).

Partial_dens depends on y(11).

Thus dydz(13) should depend at least on y(10), y(11), dydz(10),dydz(11) and molar_fractions from which I don't know how they are computed from y(1),...,y(12).

Is it not possible for you to write down dens_mix as a function of y(1),...,y(12) and to differentiate this equation with respect to z ?

If you have difficulties, make it symbolically with MATLAB.

[Gal Shaked]

the molar fraction is as follows:

%%molar fractions (No units)
molar_fractions = [0 0 0 0];
for c = 1:length(molar_fractions)
    molar_fractions(c) = y(c*2)/c_total;
end

that means, it depends on y(2), y(4), y(6), y(8)- those are solved with second order ODE's as you can see equations 2 to 9.

i will try to write down dens_mix as a function of y(1),...,y(12), but did you mean to differentiate the equation manually?

[Torsten]

As said: If you can't do it, use MATLAB's symbolic toolbox.

Everything would be much easier if you used an initial-value instead of a boundary value solver, but - as you said - this is not the case unfortunately.

[Gal Shaked]

i understand. that's became too complicated:

dydz(13)=(dydz(1))*(y(10)/(Rconstant*y(11)*y(1)))*(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))));

and still it's constant

how complicated is it to change the code and use symbolic toolbox in your opinion? is it means that i need to change the code comletely?

[Torsten]

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]

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.

Answer 2

[Torsten]

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) = 

[Torsten]

I must admit that I don't know why you need a differential equation for u = y(1).

From

d(u*rho) = 0 

you get

u*rho = constant 

for all z, thus

u(z) = (u*rho)(@z=0) / rho(z) = (m_dot_in/A_in) / rho(z)

[Gal Shaked]

thanks again

first of all i copy-pasted you previous comment of ODE's into the code and it couldnt run due to jacobian error.

in the picture below i showed how i use overall mass balance.

i didn't fully understand how you recommend solving the velocity u without a differential equation, and im not sure it's going to help because i need to keep the equation in the code as general as possible. that why i prefer working with that equation.

the thing is, when i write a display line somewhere in the code of the dens_mix :

disp(dens_mix);

the code displays all the densities it calculates during the numerical analysis. im not sure but can i use it somehow as a function in order to use it with gradient or diff ? i mean maybe without adding a new ODE dydz(13) ? i am not sure if it's possible because the simultaneously calculation

[Torsten]

I gave you the derivative of dens_mix with respect to z. I don't understand what you need more to calculate dydz(1).

I also told you that a 13th equation is superfluous because you now know the derivative of dens_mix with respect to z.

The only thing you have to do is copy the code from above and run it in MATLAB.

But as said: According to the continuity equation, rho*u = constant over the length z (assuming that the cross section of the tube remains constant). You know the inflowing mass flow, you know density_mix at position z - thus you know velocity u at position z from the formula I gave you.

[Gal Shaked]

hi, i know you gave me the 13th equation, the thing is- the code can't run it, it getting stuck by errors.

for now it's might be okay to use rho*u = constant, but for further use of the code it will be necessary to solve it with the differential equation, that why im not a big fan of that solution, anyway- thanks i really appreciate that.

but i have 2 question for that slving method:

1. 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
2. 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)

subplot(3,2,1);
plot(sol.x, sol.y(1, :), "blue");
grid on
title("Velocity");

[Torsten]

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]

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