필터 지우기
필터 지우기

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

조회 수: 11 (최근 30일)
Yuji Zhang
Yuji Zhang 2013년 6월 14일
댓글: Felipe Padua 2021년 10월 12일
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~

채택된 답변

Walter Roberson
Walter Roberson 2013년 6월 14일
>> 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
Yuji Zhang
Yuji Zhang 2013년 6월 15일
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 2013년 6월 15일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (2개)

Azzi Abdelmalek
Azzi Abdelmalek 2013년 6월 14일
편집: Azzi Abdelmalek 2013년 6월 14일
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
Yuji Zhang
Yuji Zhang 2013년 6월 15일
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?

댓글을 달려면 로그인하십시오.


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

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Logical에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by