How to calculate first and second derivative for concentration?

Hi, I have a code that show the concentration. I want to find the first and second derivative. How to calculate that? This is my code

 clc;clear;

%input parameter
  delta=50;
  gamma=75;
  K1=10^-4;
  K2=5*10^-4;
  K3=10^-3;
  K4=5*10^-3;
  K5=10^-2;
  K6=5*10^-2;
  Ko=0.1;
  n=6;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  mu_2=10^-3;
  mu_3=10^-3;
  mu_4=10^-3;
  mu_5=10^-3;
  mu_6=10^-3;
  mu_o=10^-4;
  mu_p= 10^-5;
  
  
%input initial condition
M1(1)=10;
M2(1)=0;
M3(1)=0;
M4(1)=0;
M5(1)=0;
M6(1)=0;
O(1)=0;
P(1)=0;

%input for time
t(1)=0;
dt=0.01; %time interval
t=0:dt:360; %time span

%input empty array
T=zeros(length(t)+1,1); %empty array for t
M1=zeros(length(t)+1,1); %empty array for M1
M2=zeros(length(t)+1,1); %empty array for M2
M3=zeros(length(t)+1,1); %empty array for M3
M4=zeros(length(t)+1,1); %empty array for M3
M5=zeros(length(t)+1,1); %empty array for M3
M6=zeros(length(t)+1,1); %empty array for M3
O=zeros(length(t)+1,1);
P=zeros(length(t)+1,1);

sumter=K2*M2+K3*M3+K4*M4+K5*M5;


%first and second derivative 
df=diff(M1(j+1),t);

Array indices must be positive integers or logical values.

d2f=diff(df,t);

for j = 1:length(t)
    T(j+1)=T(j)+dt;
    M1(j+1)=M1(j)+1./(1+exp(-T(j))); 
    
    M1(j+1) = M1(j)+(dt*(delta*M1(j+1)*(1-(M1(j+1)/gamma))-2*K1*M1(j+1)*M1(j+1)-M1(j+1)*sumter(j+1))-((Oa-n)*K6*M1(j+1)*M6(j+1))-((Pa-Oa)*Ko*M1(j+1)*O(j+1))-(mu_1*M1(j+1)));
    M2(j+1) = M2(j)+(dt*(K1*M1(j)*M1(j)-K2*M1(j)*M2(j))-(mu_2*M2(j+1)));
    M3(j+1) = M3(j)+(dt*(K2*M1(j)*M2(j)-K3*M1(j)*M3(j))-mu_3*M3(j));
    M4(j+1) = M4(j)+(dt*(K3*M1(j)*M3(j)-K4*M1(j)*M4(j))-mu_4*M4(j));
    M5(j+1) = M5(j)+(dt*(K4*M1(j)*M4(j)-K5*M1(j)*M5(j))-mu_5*M5(j));
    M6(j+1) = M6(j)+(dt*(K5*M1(j)*M5(j)-K6*M1(j)*M6(j))-mu_6*M6(j));
    O(j+1)  = O(j)+(dt*(K6*M1(j)*M6(j)-Ko*M1(j)*O(j)-mu_o*O(j)));
    P(j+1)  = P(j)+(dt*(Ko*M1(j)*O(j)-mu_p*P(j)));
end

subplot (2,2,1)
%figure
%hold on
plot(T,M1,'r','Linewidth',2)
xlabel('Time (days)')
ylabel('Concentration (gr/ml)')
title ('M1')

subplot (2,2,2)
%figure
%hold on
plot(T,M2,'g','Linewidth',2)
xlabel('Time (days)')
ylabel('Concentration (gr/ml)')
title ('M2')

subplot (2,2,3)
%figure
%hold on
plot(T,M3,'b','Linewidth',2)
xlabel('Time (days)')
ylabel('Concentration (gr/ml)')
title ('M3')

subplot (2,2,4)
%figure
%hold on
plot(T,M4,'m','Linewidth',2)
xlabel('Time (days)')
ylabel('Concentration (gr/ml)')
title ('M4')

subplot (2,2,1)
%figure
%hold on
plot(T,M5,'r','Linewidth',2)
xlabel('Time (days)')
ylabel('Concentration (gr/ml)')
title ('M5')

subplot (2,2,2)
%hold on
plot(T,M6,'g','Linewidth',2)
xlabel('Time (days)')
ylabel('Concentration (gr/ml)')
title ('M6')

subplot (2,2,3)
%figure
%hold on
plot(T,O,'b','Linewidth',2)
xlabel('Time (days)')
ylabel('Concentration (gr/ml)')
title ('O')

subplot (2,2,4)
%figure
%hold on
plot(T,P,'m','Linewidth',3)
xlabel('Time (days)')
ylabel('Concentration (gr/ml)')
title ('P')

Answer 1

[Walter Roberson]

%first and second derivative 
df=diff(M1(j+1),t);

You have that line before your for j loop, so at that point j is its default value, sqrt(-1)

M1=zeros(length(t)+1,1); %empty array for M1

M1 is a numeric array

df=diff(M1(j+1),t);

M1(j+1) [with valid scalar j index] would be a scalar numeric value. diff() applied to numeric values is numeric differences. diff(X,N) is X(2:end)-X(1:end-1), repeated N times -- which requires that the second parameter is a positive integer numeric scalar. But...

t=0:dt:360; %time span

it is non-scalar, non-positive, and not integer.

There are two major functions named diff . The one that applies to most data types is what I mentioned above, X(2:end)-X(1:end-1) . There is also diff symbolic derivative. Symbolic derivative can work in two modes:

• if the first parameter is of data type sym, symfun, or symmatrix, then it is an expression whose derivative is to be taken, and the remaining parameters describe which variable to take the derivative with respect to, and how many derivatives to take;
• if the first parameter to the symbolic diff function is numeric and the second parameter is symbolic variable, then the numeric variable is converted to symbolic equivalent and the derivative is taken with respect to the variable... but the symbolic equivalent of a numeric variable is of course constant, and the derivative of a constant with respect to a variable is 0. So in this case, diff(Numeric_expression, symbolic_variable) would end up returning zeros the same size as the Numeric_expression

You, though, are trying to use diff() with numeric first parameter and numeric second parameter, so you are not invoking the symbolic derivative function at all.

What should you do? Well, you should perhaps create your M* matrixes over time (recording the matrix at each time step) and then call gradient to estimate the numeric gradient.

But it depends... do you want to take the derivative with respect to time, like your code suggests... or do you want to take the derivative with respect to space?

Answer 2

[Torsten]

It seems the concentrations are solutions of differential equations, thus something like dy/dt = f(y,t).

So you can see that after you have solved the differential equations, you get the first derivative simply as dy/dt = f(y,t).

For the second derivative, you get

d^2y/dt^2 = d/dt (f(y,t)) = df/dy * dy/dt + df/dt = df/dy * f(t,y) + df/dt.

Further, there seems to be an error in your equations:

M2(j+1) = M2(j)+(dt*(K1*M1(j)*M1(j)-K2*M1(j)*M2(j))-(mu_2*M2(j+1)));

M2(j+1) appears on both sides of your equation. Since you initialized M2 to 0, the last term will always be 0 and so the equation is the same as if you had written

M2(j+1) = M2(j)+(dt*(K1*M1(j)*M1(j)-K2*M1(j)*M2(j)));

Other incorrect settings are the use of O(j+1) and sumter(j+1) in the equation for M1(j+1).