time derivative of pressure

I have a .txt file with two columns of values: the first is the pressure, the second is the time (at each instant of time I have an assigned value of pressure). I need to calculate the time derivative of the pressure.

P.s: the method like: dpdt = diff(p) / diff(t), with p e t the values of the .txt file, cannot be used since bring to wrong results.

The suggestion was to use backward / forward / central differences, but I do not know how to use.

Thanks!

[Torsten]

If diff(p)/diff(t) is not satisfactory, then maybe there is too much noise in the pressure measurements to be differentiated directly.

Could you supply the data ?

[Simone Pavesi]

Data are not measured, but are from a simulation.

these are a series of data: first column time, second column pressure (I know that much of the pressure values are equal one to each other, and probably this is the problem, but simulation give me 40000 values of pressure and time, for 4ms of simulation):

4.59e-05 5.99978e+06

4.6e-05 5.99977e+06

4.61e-05 5.99977e+06

4.62e-05 5.99977e+06

4.63e-05 5.99977e+06

4.64e-05 5.99977e+06

4.65e-05 5.99977e+06

4.66e-05 5.99977e+06

4.67e-05 5.99977e+06

4.68e-05 5.99977e+06

4.69e-05 5.99977e+06

4.7e-05 5.99977e+06

4.71e-05 5.99977e+06

4.72e-05 5.99977e+06

4.73e-05 5.99977e+06

4.74e-05 5.99977e+06

4.75e-05 5.99977e+06

4.76e-05 5.99977e+06

4.77e-05 5.99977e+06

4.78e-05 5.99977e+06

4.79e-05 5.99976e+06

[Star Strider]

Assuming the matrix of those values is ‘A’,

dpdt = gradient(A(:,2)) ./ gradient(A(:,1))

See the documentation on the gradient function for details.

.

[Torsten]

@Simone Pavesi

Do you know what the gradient of a constant function is ? This will be the result in your case (except for the first and last element).

Maybe outputting pressure with more digits will change the situation.

[Prakash S R]

@Simone Pavesi

What does "wrong results" mean when you do 1st difference? perhaps you are doing diff(p) / diff(t) instead of elementwise: diff(p) ./ diff(t)

If by "wrong" you mean you get NaNs, then it may be that the time column has repeated values, giving a denominator of 0. Eliminating the repeats will solve the problem. But I would examine whay they are there in the first place.

gradient will simply compute the 2nd difference (except at start/end). It will work when there is a single repeat of the t values, but if you have two or more repeats, that will return NaNs as well.

@Torsten

The gradient of a constant should legitimately return 0. The problem here may be that the time column has repeated values, resulting in NaN when you do gradient

[Simone Pavesi]

@Torsten I know what is the output of a constant function gradient. The problem is not the single derivative result, but how to write it, I mean, how to write it with the backward/forward/central difference

[Simone Pavesi]

@Star Strider thanks for your suggestion, unfortunately in the way you seggest does not work.

[Torsten]

Forward scheme: (p(2:end)-p(1:end-1))./(t(2:end)-t(1:end-1)) (gives derivative at t(1),...t(end-1))

Central scheme: (p(3:end)-p(1:end-2))./(t(3:end)-t(1:end-2)) (gives derivative at t(2),...,t(end-1))

Backward scheme: (p(2:end)-p(1:end-1))./(t(2:end)-t(1:end-1)) (gives derivative at t(2),...t(end))

[Star Strider]

@Simone Pavesi — It works. The data need more precision than is posted —

A = [4.59e-05	5.99978e+06
    4.6e-05	5.99977e+06
    4.61e-05	5.99977e+06
    4.62e-05	5.99977e+06
    4.63e-05	5.99977e+06
    4.64e-05	5.99977e+06
    4.65e-05	5.99977e+06
    4.66e-05	5.99977e+06
    4.67e-05	5.99977e+06
    4.68e-05	5.99977e+06
    4.69e-05	5.99977e+06
    4.7e-05	    5.99977e+06
    4.71e-05	5.99977e+06
    4.72e-05	5.99977e+06
    4.73e-05	5.99977e+06
    4.74e-05	5.99977e+06
    4.75e-05	5.99977e+06
    4.76e-05	5.99977e+06
    4.77e-05	5.99977e+06
    4.78e-05	5.99977e+06
    4.79e-05	5.99976e+06];

dpdt = gradient(A(:,2)) ./ gradient(A(:,1));

figure
yyaxis left
plot(A(:,1), A(:,2))
ylabel('$P(t)$', 'Interpreter','latex')
yyaxis right
plot(A(:,1), dpdt)
ylabel('$\frac{dP(t)}{dt}$', 'Interpreter','latex')
xlabel('$t$', 'Interpreter','latex')

.

[Simone Pavesi]

Thank you at all guys!