Calculating derivative of discrete data

I have x,y data which describes change in diameter (y) over time (x) (time, diameter). I would like to find the derivative dD/dt in order to obtain the viscosity (Eq.1), and then plot against strain (Eq. 2). I am testing it out with some data that I grabbed from a paper, however, I cannot seem to replicate their results. Sigma is a constant.

Eq. 1

Eq.2

I am unsure if I need to do a spline fit of the data or whether there is an issue with me calculating the derivative. Either way there is some "noise" in my results. Any help would be great!

My basic code is below:

clearvars
close all

t = [0.0283,	0.1414,	0.3563,	0.4808,	0.6052,	0.7296,	0.8428,	0.9615,	1.0803,	1.1991,	1.3292,	1.4536, 1.5667,	1.6799,	1.8043,	1.9061,	2.0249,	2.2851,	2.5396,	2.7376,	2.8903,	3.0090,	3.3201,	3.5068,	3.7387,	4.0498,	4.2081];
D =[0.2553,	0.2287,	0.1793,	0.1422,	0.0960,	0.0618,	0.0541,	0.0508,	0.0474,	0.0449,	0.0427,	0.0405,0.0382,	0.0371,	0.0348,	0.0340,	0.0310,	0.0289,	0.0268,	0.0239,	0.0215,	0.0204,	0.0177,	0.0154,	0.0138,	0.0125,	0.0112];

extenData_strain = [3.05,	3.34,	3.49,	3.65,	3.83,	4.03,	4.20,	4.38,	4.56,	4.74,	4.92,	5.10,	5.29,	5.47,	5.66,	5.86,	6.10,	6.36];
extenData_vis = [1.03,	2.65,	3.27,	3.70,	4.20,	4.65,	5.04,	5.45,	5.83,	6.16,	6.52,	6.80,	7.11,	7.39,	7.58,	7.75,	7.82,	7.84];


N = length(t);
strain = 2*log(D(1) ./ D); %calcualte hencky strain
dD_D = diff(D); %differentiate diameter 
dt = diff(t);
der = dD_D./dt;
der(N) = der(N-1);

sigma=0.066;
extensionalVis = -sigma./der;

figure(1)
semilogy(t, D, 'o')
hold on 
xlabel("Time(ms)")
ylabel("Diameter (mm)")
xlim([0 5])
ylim([0.01 0.3])
set(gca, 'YScale', 'log')

figure(2)
plot(strain, extensionalVis, 'o', 'DisplayName','My results')
hold on 
plot(extenData_strain,extenData_vis, 'x', 'DisplayName', 'Paper results' )
xlabel("Hencky strain")
ylabel("Extensional viscosity (Pa.s)")
xlim([3 6.5])
ylim([0 9])
legend();

[Mathieu NOE]

hi

well, it's BIG noise , if we still consider that as noise

tryed a few options with interpolation then smoothdata and exponential fit , maybe there are other options but your data are quite off from the theory so I guess no tool will give you the expected results unless you apply a massive amount of corretion , up to the point we could consider that as "data manipulation"

clearvars
close all
t = [0.0283,	0.1414,	0.3563,	0.4808,	0.6052,	0.7296,	0.8428,	0.9615,	1.0803,	1.1991,	1.3292,	1.4536, 1.5667,	1.6799,	1.8043,	1.9061,	2.0249,	2.2851,	2.5396,	2.7376,	2.8903,	3.0090,	3.3201,	3.5068,	3.7387,	4.0498,	4.2081];
D =[0.2553,	0.2287,	0.1793,	0.1422,	0.0960,	0.0618,	0.0541,	0.0508,	0.0474,	0.0449,	0.0427,	0.0405,0.0382,	0.0371,	0.0348,	0.0340,	0.0310,	0.0289,	0.0268,	0.0239,	0.0215,	0.0204,	0.0177,	0.0154,	0.0138,	0.0125,	0.0112];

tt = linspace(min(t),max(t),150);
DD = interp1(t,D,tt);
t = tt;
D = DD;

extenData_strain = [3.05,	3.34,	3.49,	3.65,	3.83,	4.03,	4.20,	4.38,	4.56,	4.74,	4.92,	5.10,	5.29,	5.47,	5.66,	5.86,	6.10,	6.36];
extenData_vis = [1.03,	2.65,	3.27,	3.70,	4.20,	4.65,	5.04,	5.45,	5.83,	6.16,	6.52,	6.80,	7.11,	7.39,	7.58,	7.75,	7.82,	7.84];
N = length(t);
strain = 2*log(D(1) ./ D); %calcualte hencky strain
% dD_D = diff(D); %differentiate diameter 
dt = mean(diff(t));
% der = dD_D./dt;
% der(N) = der(N-1);
der = gradient(D)/dt;
sigma=0.066;
extensionalVis = -sigma./der;


% crude exponential fit to data
    ind = find(strain>=3);
    extensionalVis = extensionalVis(ind);
    strain = strain(ind);

    me = extensionalVis(end);
%     me = max(extensionalVis);
    tmp = me - extensionalVis;
    threshold = 1e-1;
    tmp(tmp<threshold) = threshold;

    P = polyfit(strain, log(tmp), 1);
    m = P(1);
    b = exp(P(2));
    tmp2 = b*exp((strain)*m);
    extensionalVis2 = me - tmp2;
    
% smooth the data
extensionalVis3 = smoothdata(extensionalVis,'lowess',75);

figure(1)
semilogy(t, D, 'o')
hold on 
xlabel("Time(ms)")
ylabel("Diameter (mm)")
xlim([0 5])
ylim([0.01 0.3])
% set(gca, 'YScale', 'log')
figure(2)
plot(strain, extensionalVis,'o', strain, extensionalVis2,'.-', strain, extensionalVis3,'.-')
legend({'My results';'exp fit';'smoothed'});hold on 
plot(extenData_strain,extenData_vis, 'x', 'DisplayName', 'Paper results' )
xlabel("Hencky strain")
ylabel("Extensional viscosity (Pa.s)")
% xlim([3 6.5])
% ylim([0 9])
legend();

[Mathieu NOE]

Another try ....

clearvars
close all
t = [0.0283,	0.1414,	0.3563,	0.4808,	0.6052,	0.7296,	0.8428,	0.9615,	1.0803,	1.1991,	1.3292,	1.4536, 1.5667,	1.6799,	1.8043,	1.9061,	2.0249,	2.2851,	2.5396,	2.7376,	2.8903,	3.0090,	3.3201,	3.5068,	3.7387,	4.0498,	4.2081];
D =[0.2553,	0.2287,	0.1793,	0.1422,	0.0960,	0.0618,	0.0541,	0.0508,	0.0474,	0.0449,	0.0427,	0.0405,0.0382,	0.0371,	0.0348,	0.0340,	0.0310,	0.0289,	0.0268,	0.0239,	0.0215,	0.0204,	0.0177,	0.0154,	0.0138,	0.0125,	0.0112];
extenData_strain = [3.05,	3.34,	3.49,	3.65,	3.83,	4.03,	4.20,	4.38,	4.56,	4.74,	4.92,	5.10,	5.29,	5.47,	5.66,	5.86,	6.10,	6.36];
extenData_vis = [1.03,	2.65,	3.27,	3.70,	4.20,	4.65,	5.04,	5.45,	5.83,	6.16,	6.52,	6.80,	7.11,	7.39,	7.58,	7.75,	7.82,	7.84];
N = length(t);
strain = 2*log(D(1) ./ D); %calcualte hencky strain
dD_D = diff(D); %differentiate diameter 
dt = diff(t);
der = dD_D./dt;
der(N) = der(N-1);
sigma=0.066;
extensionalVis = -sigma./der;

    % another try 
    strain(extensionalVis>10) = []; % remove this large outlier
    extensionalVis(extensionalVis>10) = [];
    ind = find(strain>=3);
    extensionalVis = extensionalVis(ind);
    strain = strain(ind);
    
    strain2 = linspace(min(strain),max(strain),100);
    extensionalVis2 = interp1(strain,extensionalVis,strain2);
    strain = strain2;
    extensionalVis = extensionalVis2;

% crude exponential fit to data
    me = extensionalVis(end);
%     me = max(extensionalVis);
    tmp = me - extensionalVis;
    threshold = 3e-1;
    tmp(tmp<threshold) = threshold;

    P = polyfit(strain, log(tmp), 1);
    m = P(1);
    b = exp(P(2));
    tmp2 = b*exp((strain)*m);
    extensionalVis2 = me - tmp2;
    
% smooth the data
extensionalVis3 = smoothdata(extensionalVis,'lowess',75);

figure(1)
semilogy(t, D, 'o')
hold on 
xlabel("Time(ms)")
ylabel("Diameter (mm)")
xlim([0 5])
ylim([0.01 0.3])
% set(gca, 'YScale', 'log')
figure(2)
plot(strain, extensionalVis,'o', strain, extensionalVis2,'.-', strain, extensionalVis3,'.-')
legend({'My results';'exp fit';'smoothed'});hold on 
plot(extenData_strain,extenData_vis, 'x', 'DisplayName', 'Paper results' )
xlabel("Hencky strain")
ylabel("Extensional viscosity (Pa.s)")
% xlim([3 6.5])
% ylim([0 9])
legend();

[N_mat]

Just to clarify. The 'D' and 't' values are also from the paper. The authors calculate the extensional viscosity from 'D' and 't'. So either there is something off about the data or it is my method. I am using their data to validate my method.

[Mathieu NOE]

ok

something is wrong somewhere

can you share the publication ?

[N_mat]

Apologies, I have just read a section in the their methods. "In order to calculate the extensional viscosity, the diameter decay was fit with a spline and then differentiated as described below". It sounds that they differentiated the spline, this might reduce the noise?

I have attached the paper. I have been looking at the diameter data in Figure 2a (circle data) and comparing with the extensional viscosity data in Figure 4a (circle data)

[Mathieu NOE]

Just finished a first quick reading but seems the fit process is not much detailled

I see the D vs t curve must be splitted in two regions and then fitted separately ?

[Mathieu NOE]

still I feel their fitted model (for D vs t) are quite close to the actual measurements - so why the code generate so much outliers (with raw measurements) ... is still strange to me.

[Mathieu NOE]

tried some spline fitting (with this FEX : SPLINEFIT - File Exchange - MATLAB Central (mathworks.com) )

but honestly, could not really make the final plot much better (tested with n = 6 and 11 pieces for the spline fit)

clearvars
close all
t = [0.0283,	0.1414,	0.3563,	0.4808,	0.6052,	0.7296,	0.8428,	0.9615,	1.0803,	1.1991,	1.3292,	1.4536, 1.5667,	1.6799,	1.8043,	1.9061,	2.0249,	2.2851,	2.5396,	2.7376,	2.8903,	3.0090,	3.3201,	3.5068,	3.7387,	4.0498,	4.2081];
D =[0.2553,	0.2287,	0.1793,	0.1422,	0.0960,	0.0618,	0.0541,	0.0508,	0.0474,	0.0449,	0.0427,	0.0405,0.0382,	0.0371,	0.0348,	0.0340,	0.0310,	0.0289,	0.0268,	0.0239,	0.0215,	0.0204,	0.0177,	0.0154,	0.0138,	0.0125,	0.0112];
extenData_strain = [3.05,	3.34,	3.49,	3.65,	3.83,	4.03,	4.20,	4.38,	4.56,	4.74,	4.92,	5.10,	5.29,	5.47,	5.66,	5.86,	6.10,	6.36];
extenData_vis = [1.03,	2.65,	3.27,	3.70,	4.20,	4.65,	5.04,	5.45,	5.83,	6.16,	6.52,	6.80,	7.11,	7.39,	7.58,	7.75,	7.82,	7.84];
N = length(t);


% Spline fit D vs t : 
% Breaks interpolated from data
nb = 6; % 6 or 11 ? 
pp2 = splinefit(t, D,nb);  % nb + 1 breaks, nb pieces
D2 = ppval(pp2,t);

figure(1)
semilogy(t, D, 'o',t, D2)
hold on 
xlabel("Time(ms)")
ylabel("Diameter (mm)")
xlim([0 5])
ylim([0.01 0.3])
% set(gca, 'YScale', 'log')

%% next
D = D2;
strain = 2*log(D(1) ./ D); %calcualte hencky strain
dD_D = diff(D); %differentiate diameter 
dt = diff(t);
der = dD_D./dt;
der(N) = der(N-1);


sigma=0.066;
extensionalVis = -sigma./der;

figure(2)
plot(strain, extensionalVis, 'o', 'DisplayName','My results')
hold on 
plot(extenData_strain,extenData_vis, 'x', 'DisplayName', 'Paper results' )
xlabel("Hencky strain")
ylabel("Extensional viscosity (Pa.s)")
xlim([3 6.5])
ylim([0 9])
legend();