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

25 Ansichten (letzte 30 Tage)
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~

Akzeptierte Antwort

Walter Roberson
Walter Roberson am 14 Jun. 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 Kommentare
Yuji Zhang
Yuji Zhang am 15 Jun. 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 am 15 Jun. 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.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (2)

Azzi Abdelmalek
Azzi Abdelmalek am 14 Jun. 2013
Bearbeitet: Azzi Abdelmalek am 14 Jun. 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 Kommentar
Yuji Zhang
Yuji Zhang am 15 Jun. 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?

Melden Sie sich an, um zu kommentieren.


Iain
Iain am 14 Jun. 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 Kommentare
Iain
Iain am 15 Jun. 2013
Two consecutive points on your curve "y" :P
Maybe I should have them put my name in all caps...
Felipe Padua
Felipe Padua am 12 Okt. 2021
You could also use
middle_of_the_time_between_y1_and_y2 = movemean(t, 2, 'endpoints', 'discard')

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

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

Start Hunting!

Translated by