Numerical differentiation in time in function used by ODE

Hi all, I'm working on a code that uses an ODE integration instead of a Runge - Kutta method (4th order).

By now I'm trying to see if it's possibile to numerically differentiate in time a variable calculated inside a function that is then given as an input to ODE. Essentially my function will compute values at single instants as time proceeds and for this reason I can't define an array on which I can compute a numerical differentiation with diff function.

I'm attaching an image of the full equation:

The variable I want to differentiate numerically is .

Thanks in advance.

[Torsten]

Which ODE integrator do you want to use ? Is rho_l_dot a solution variable ?

[Mithrandir]

@Torsten ODE45 or ODE113. is not inserted as a solution of the integration but computed through some variables that results from integration.

[James Tursa]

Show us how rho_dot_L is computed.

[Mithrandir]

@James Tursa By now I'm using chain rule to compute this derivative according to this formula:

considering that density has a dependency on pressure p and temperature T.

In my model I compute derivative of p and T in time according to conservation of energy equation and definition of density.

[Torsten]

Do you want to differentiate rho_l or rho_l_dot with respect to time ? You wrote you want to differentiate rho_l_dot again to get rho_l_doubledot. But maybe that was wrong ?

[Mithrandir]

@Torsten Only rho_l since no second derivative will be needed. Sorry for stating it wrong.

Answer 1

[Torsten]

ode15s can solve systems of differential equations of the form M*y' = f(t,y) where M is the so called mass matrix. If a row j of this matrix is zero, the corresponding equation reads 0 = f_j(t,y) , thus is an algebraic equation.

In your case, define a new algebraic solution variable v as v = rho_l and add the algebraic equation v - rho_l = 0 to your system of differential-algebraic equations to be solved. If you now have to use dv/dt (= drho_l/dt) somewhere in your code, you can include it in your computations by setting a coefficient different from 0 in the mass matrix M of the system in the column that corresponds to v .

[Mithrandir]

@Torsten Thanks Torsten, I'll try this way. I'm new to Matlab but thanks a lot for the advice.

[Torsten]

I adjusted the answer for that you want rho_l_dot, not rho_l_doubledot.

And you should start with a simple test example to see how it works to implicitly get the time derivative of an algebraic variable in the way I described.

[Mithrandir]

@Torsten Thanks Torsten, much appreciated.