Numerical derivative in matlab

Hello,

I'm using Matlab 2012 and would need the numerical derivative of a function. What is the Malab command for that ?

Thank you in advance

[Leah]

I would start here http://en.wikipedia.org/wiki/Numerical_differentiation

I like your user name

Answer 1

[Azzi Abdelmalek]

If you have a series of data

t=[0 1 2 3 4 5] % Vector time
y=[1 2 4 8 12 13] % Your signal
% The derivative can be approximated :

diff(y)./diff(t)

%You can also use symbolic derivative

syms t y
y=cos(t)
dy=diff(y)

[dont panic]

thank you for your reply :) !

Answer 2

[Rohit Garud]

You can use the gradient() function in newer versions of matlab

dy_by_dt = gradient(y(:)) ./ gradient(t(:))

[Daphne PARLIARI]

Indeed this works like a charm.

Do you have any idea how to convert dy_by_dt to percentage change?

[Walter Roberson]

percentage_change = dy_by_dt ./ y(:) * 100;

[Daphne PARLIARI]

I have the attached data, which are daily values of temperature and deaths. I want to calculate percent increase in mortality per 1oC increase above a certain heat threshold. Eg. if the threshold is 30oC, there is a 2.1% increase in mortality per 1oC increase in temperature (31oC, 32oC, ...).

Currently I have written this code although I am not sure it works the way it should:

idx = Daily_Data.Daily_T(:) > 30;
Daily_Data_new  = Daily_Data(idx,:);
dy_by_dt = gradient(Daily_Data_new.Daily_Deaths(:)) ./ gradient(Daily_Data_new.Daily_T(:));
percentage_change = dy_by_dt ./ Daily_Data_new.Daily_Deaths(:) * 100;

I would appreciate any help...

[Walter Roberson]

I don't think that is right at all.

Daily_Data = readtable('Daily_Data_MoT.xlsx');

idx = Daily_Data.Daily_T(:) > 30;
Daily_Data_new  = Daily_Data(idx,:);
[G, Temperature] = discretize(Daily_Data_new.Daily_T, 30:37);
d_count = accumarray(G, Daily_Data_new.Daily_Deaths, [length(Temperature), 1]);
d_mean = accumarray(G, Daily_Data_new.Daily_Deaths, [length(Temperature), 1], @mean);

bar(Temperature, d_count); xlabel('temperature (C)'); ylabel('Total deaths')

bar(Temperature, d_mean); xlabel('temperature (C)'); ylabel('average deaths')

You could do some kind of statistical analysis to show that the mean deaths increases with temperature 30 to 32, but it decreases after that up to 34, and then increases again at 36 but not to as much as 32.

[Walter Roberson]

You could look at gradient(d_mean) but I do not think you will find it enlightening.

[Daphne PARLIARI]

Wow, the graphs were great, thank you!

Could you please explain what does every term in

d_count = accumarray(G, Daily_Data_new.Daily_Deaths, [length(Temperature), 1]);

do and what is the difference with

d_mean = accumarray(G, Daily_Data_new.Daily_Deaths, [length(Temperature), 1], @mean) ?

In particular, I don't quite get what G is doing.

[Walter Roberson]

G is a grouping number, similar in spirit to using findgroups(), except that findgroups() treats every unique value as a different bin but you want to treat ranges as different bins.

[G, Temperature] = discretize(Daily_Data_new.Daily_T, 30:37);

is pretty much the same as

G = floor(Daily_Data_new.Daily_T - 30) + 1;

in that it maps a degree at a time, [30, 31) -> 1, [31,32) -> 2, [32,33) -> 3 and so on. It is the same thing exactly as histogram bin number doing a histogram of temperatures above 30 one degree per bin.

d_count = accumarray(G, Daily_Data_new.Daily_Deaths, [length(Temperature), 1]);

That code is the same as

d_count = zeros(length(Temperature),1);
DD = Daily_Data_new.Daily_Deaths;
for K = 1 : length(G)
    g = G(K);
    d_count(g) = d_count(g) + DD(K);
end

That is, it uses the G entries as index numbers into the array, and accumulates all of the Daily_Death values corresponding to the same index. Another way of looking at it would be

d_count = zeros(length(Temperature),1);
for g = 1 : 7
    d_count(g) = sum(DD(G == g));
end

and

d_mean = accumarray(G, Daily_Data_new.Daily_Deaths, [length(Temperature), 1], @mean)

would correspond to

d_mean = zeros(length(Temperature),1);
for g = 1 : 7
    d_mean(g) = mean(DD(G == g));
end

and in turn you could rewrite the code to:

Temperature = 30:37;
nt = length(Temperature)-1;
d_count = zeros(nt,1);
d_mean = zeros(nt,1);
DD = Daily_Data_new.Daily_Deaths;
DT = Daily_Data_new.Daily_T;
for g = 1 : nt
    d = DD(DT >= Temperature(g) & DT < Temperature(g+1));
    d_count(g) = sum(d);
    d_mean(g) = mean(d);
end

[Daphne PARLIARI]

Thanks for the detailed answer, that made all clear! Let me pose another question, if you have the time of course.

I want to calculate percent change with reference to a constant temperature, e.g. 30oC, and add that percentage on top of each column. When I try

x = b.XData;
y = b.YData;
reference_temp = y(1);  % 30oC in plot
percent_change = 100 * (y - reference_temp) / reference_temp;
for k = 1:length(percent_change)
    text(x(k+1),y(k+1),sprintf('%0.1f%%',percent_change(k)),...
        'HorizontalAlignment','center',...
        'VerticalAlignment','bottom')
end

I am getting the Error: Index exceeds the number of array elements (8).

[Walter Roberson]

Is b the handle of the bar() ?

length(percent_change) is going to be the same as length(y) . You loop k 1 to length percent_change which is the same as loop k 1 to length y. You access y(k+1) but when k becomes length(y) then you are accessing y(length(y)+1) which is after the end of y.

[Daphne PARLIARI]

Oh I got confused... What would be the solution to that?

[Walter Roberson]

for k = 2:length(percent_change)
  text(x(k),y(k),sprintf('%0.1f%%',percent_change(k)),...
      'HorizontalAlignment','center',...
      'VerticalAlignment','bottom')
end

Answer 3

[Tamas Kis]

https://www.mathworks.com/matlabcentral/fileexchange/89719-numerical-differentiation-differentiate

Function that can differentiate a set of points given as vectors (such as how you would use cumtrapz to integrate), or differentiate a function defined using a function handle. It can differentiate at any specified point(s) or over an entire domain. Also is of slightly higher accuracy than using diff(y)./diff(t) because it uses central approximation at the interior nodes.