Effacer les filtres
Effacer les filtres

doing derivative using diff(Y)/dT makes the vector shorter

11 vues (au cours des 30 derniers jours)
Yuji Zhang
Yuji Zhang le 14 Juin 2013
Commenté : Felipe Padua le 12 Oct 2021
Hi everybody,
I'm doing derivative of a curve Y-T I think it's:
T = linspace(-t, t, n); Y = somefuction; dT = T(2)-T(1); DY1 = diff(Y)/dT;
But then DY1 is one element shorter than Y. How do people usually deal with this?
I'm currently dealing with this by shorten the T axis:
plot( T(1:end-1), DY1 );
I don't know whether this is the best way... Is there are relative standard way? Let me know. Thanks everyone~

Réponse acceptée

Walter Roberson
Walter Roberson le 14 Juin 2013
>> dT = 1;
>> x = 1:dT:3;
>> y = x.^2;
>> diff(y) ./ dT
3 5
>> diff('x^2', 'x')
2 * x
>> subs( diff('x^2', 'x'), 'x', x(1:length(y)) )
2 4
Thus we can see that using diff(y)/dT does not give us the same results as if we worked symbolically.
Question then: at what x values are [3 5] the correct derivative according to symbolic methods ? 2 * xp = [3 5]... by examination, xp must be [3/2 5/2]. Which, by no coincidence at all is the evaluation at the midpoints between x(n) and x(n+1) -- (x(n) + x(n+1))/2, or compactly (x(1:end-1) + x(2:end))/2
If we step back for a few seconds, we can see that using the numeric formula diff(y) ./ dT assigns the entire difference y(n) to y(n+1) as if it were at x(n), but that is not how derivatives work: derivatives are the tangent around x(n) and so y(n-1) must be taken into account, not just y(n) and y(n+1). Easiest resolution is to use x(n) and x(n+1) and y(n) and y(n+1) to construct the slope associated with the midpoint (x(n) + x(n+1))/2
plot( (T(1:end-1)+T(2:end))/2, DY1 )
If, however, you need to a slope at each x(n), then you have problems with the definition of slope at the endpoints. You might, in that case, wish to use the definition predefined:
plot( T, gradient(T) )
  2 commentaires
Yuji Zhang
Yuji Zhang le 15 Juin 2013
Hi Walter~
Thanks for sharing!
I see, for 1D curve Y-T, this
>>plot( (T(1:end-1)+T(2:end))/2, DY1 );
reflects the definition of derivative strictly.
--------
Actually gradient(Y) is a vector with same length as T - by the math definition of gradient. So
>>plot ( T, gradient(Y) );
should be good. The problem of this approach is, the 1st and last element of gradient(Y) are of error. - which makes sense according to its math definition too.
For my specific purpose, I deal with a long enough numerical curve (>1000 element, not symbolic function). So I think both approach works. The second might be more convenient of the boundary values are not important.
What do you think Walter? Any discussion is appreciated everyone~
Walter Roberson
Walter Roberson le 15 Juin 2013
I have no recommendation. Both approaches are valid in different situations. When a task requires a derivative at every point, I study why it requires that in the circumstances, and use whatever endpoint calculations are most suitable for the circumstances. More often, perhaps, I would use the interior points only, in order to avoid the problem.

Connectez-vous pour commenter.

Plus de réponses (2)

Azzi Abdelmalek
Azzi Abdelmalek le 14 Juin 2013
Modifié(e) : Azzi Abdelmalek le 14 Juin 2013
Edit
(x(2)-x(1))/(t(2)-t(1)) correspond to the approximative right derivative at the point(t(1),x(1)). The last point is (t(n-1),x(n-1)), which means that you are doing right
  1 commentaire
Yuji Zhang
Yuji Zhang le 15 Juin 2013
Hi Azzi~
Yes, it's just inconvenient cause you need to worry about the length... I think gradient(Y) could be better. What do you think?

Connectez-vous pour commenter.


Iain
Iain le 14 Juin 2013
How I normally do it:
average_slope_between_y1_and_y2 = diff(Y)./diff(t);
middle_of_the_time_between_y1_and_y2 = (t(2:end)+t(1:end-1))./2;
Alternatively, fit a curve to your data, and differentiate that.
  3 commentaires
Iain
Iain le 15 Juin 2013
Two consecutive points on your curve "y" :P
Maybe I should have them put my name in all caps...
Felipe Padua
Felipe Padua le 12 Oct 2021
You could also use
middle_of_the_time_between_y1_and_y2 = movemean(t, 2, 'endpoints', 'discard')

Connectez-vous pour commenter.

Catégories

En savoir plus sur Logical dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by