how to calculate the differentiation by diff command ?

I need to calculate the derivative and second derivative of phid and thetad (w.r.t time) for the code given below. Can anyone plz help in calculating this?

m = 0.65;

d = 7.5*10^-7;

l = 0.23;

Jx = 7.5 * 10^-3;

Jy = 7.5 * 10^-3;

Jz = 1.3 * 10^-2;

b = 3.13 * 10^-5;

a1 = (Jy - Jz)/Jx ; b1 = 1/Jx;

a2 = (Jz - Jx)/Jy; b2 = 1/Jy;

a3 = (Jx - Jy)/Jz; b3 = 1/Jz;

g0 = 9.81;

c1 = 1; c3 = 1; c5 = 1; c7 = 1; c9 = 1; c11 = 1;

c2 = 1; c4 = 1; c6 = 1; c8 = 5; c10 = 1; c12 = 1;

x1(1) = 0; %% roll

x2(1) = 0;

x3(1) = 0; %% pitch

x4(1) = 0;

x5(1) = 0; %% yaw

x6(1) = 0;

x7(1) = 0; %% z position

x8(1) = 0;

x9(1) = 0; %% x poition

x10(1) = 0;

x11(1) = 0; %% y position

x12(1) = 0;

dt = 0.1;

t = 0:dt:60;

for n = 1: length(t)

phid(1) = 0;

thetad(1) = 0;

xd(:,n) = [phid(n); thetad(n); 0; 0; 0; 0; zdes; diff(zdes,t); xdes; diff(xdes,t); ydes; diff(ydes,t)];

xdd(:,n) = [0; 0; 0; 0; 0; 0; diff(zdes,t); diff(diff(zdes,t)); diff(xdes,t); diff(diff(xdes,t)); diff(ydes,t); diff(diff(ydes,t))];

xddd(:,n) = [0; 0; 0; 0; 0; 0; diff(diff(zdes,t)); diff(diff(diff(zdes,t))); diff(diff(xdes,t)); diff(diff(diff(xdes,t))); diff(diff(ydes,t)); diff(diff(diff(ydes,t)))];

e1(:,n) = phid(n) - x1(n);

e3(:,n) = thetad(n) - x3(n);

e5(:,n) = xd(5,n) - x5(n);

e7(:,n) = xd(7,n) - x7(n);

e9(:,n) = xd(9,n) - x9(n);

e11(:,n) = xd(11,n) - x11(n);

e2(:,n) = x2(n) - xdd(1,n) - c1*e1(n);

e4(:,n) = x4(n) - xdd(3,n) - c3*e3(n);

e6(:,n) = x6(n) - xdd(5,n) - c5*e5(n);

e8(:,n) = x8(n) - xdd(7,n) - c7*e7(n);

e10(:,n) = x10(n) - xdd(9,n) - c9*e9(n);

e12(:,n) = x12(n) - xdd(11,n) - c11*e11(n);

U1(n) = ( m / ( cos( x1(n) ) * cos(x3(n)) ) * (g0 + xdd(8,n) + e7(n) - c8*e8(n)));

Ux(n) = (m/(U1(n)))*(xdd(10,n) + e9(n) - c10*e10(n));

Uy(n) = (m/(U1(n)))*(xdd(12,n) + e11(n) - c12*e12(n));

phid(n+1) = asin(Ux(n)*sin(xd(5,n)) - Uy(n)*cos(xd(5,n)));

thetad(n+1) = asin( ( Ux(n)*sin(xd(5,n)) + Uy(n)*cos(xd(5,n))) )/sqrt(1-( Ux(n)*sin(xd(5,n) - Uy(n)*cos(xd(5,n))) )^2 ) ;

U2(n) = (1/b1)*(- a1*x4(n)*x6(n) + xdd(2,n) + (phid(n)-x1(n)) - c2*e2(n));

U3(n) = (1/b2)*(- a2*x2(n)*x6(n) + xdd(4,n) + (thetad(n) - x3(n)) - c4*e4(n));

U4(n) = (1/b3)*(- a3*x2(n)*x4(n) + xdd(5,n) + e5(n) - c6*e6(n));

x1(n+1) = x1(n) + dt * (x2(n));

x2(n+1) = x2(n) + dt * (a1*x4(n)*x6(n) + b1*U2(n));

x3(n+1) = x3(n) + dt * (x4(n));

x4(n+1) = x4(n) + dt * (a2*x2(n)*x6(n) + b2*U3(n));

x5(n+1) = x5(n) + dt * (x6(n));

x6(n+1) = x6(n) + dt * (a3*x2(n)*x4(n) + b3*U4(n));

x7(n+1) = x7(n) + dt * (x8(n));

x8(n+1) = x8(n) + dt * ((1/m)*(U1(n)*cos(x1(n)) * cos(x3(n))) - g0);

x9(n+1) = x9(n) + dt * (x10(n));

x10(n+1) = x10(n) + dt * ((Ux(n)* U1(n)) / m);

x11(n+1) = x11(n) + dt * (x12(n));

x12(n+1) = x12(n) + dt * ( (U1(n)*Uy(n))/m );

end

[Torsten]

Where do you use the input x vector ?

Or should x1, x2 etc. be x(1), x(2) etc. ?

And Ux(t), Uy(t) will throw an error. Functions are defined either as function handles or as symbolic expressions.

[Pallov Anand]

thanks for reply....the previous code was wrong...i just edited it now...can you plz have a look now

[Torsten]

I'm unable to make sense of your code.

Please state your problem in a mathematical notation (with formulas and without MATLAB code).

Especially I don't know what xdes, ydes and zdes are.

[Pallov Anand]

https://ieeexplore.ieee.org/abstract/document/7794007

This is the paper im trying to implement....and im again sorry in the above code regarding xd(:,n), xdd(:,n) and xddd(:,n), below is the required code for it

xd(:,n) = [phid(n); 0; thetad(n); 0; 0; 0; 0.7; 0; 2; 0; 5; 0]; % Desired position psi, z, x, y in order

xdd(:,n) = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]; % Derivative of desired position

xddd(:,n) = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]; % Double derivative of desired position

So, if you can refer to the paper,

\dot(x1),...\dot(x12), I have written as x1(n+1),...x12(n+1)

Here, I have slight modifications in my code in terms of expression for e1 to e12

xd, xdd, xddd are desired values, first derivative of the desired values and second derivative of the desired values.

In paper, phid and thetad is not given which is required.....and i have included the expressions of both in the code.

What is happening is, phid and thetad is being calculated using Ux, Uy and then xd(:,1) is taking the updated phid value.

But currently, in xdd which is first derivative of xd, I have put all zeros....but in xdd(:,1), I should derivate phid and that value should be there rather than 0.

Calculating the derivative of phid is what Im unable to do.

I hope, I can give you some idea. If not, let me know.

Thanks

[Torsten]

You just wrote down the derivatives in the problem formulation:

dphi/dt       = x2
d^2phi/dt^2   = a1*x4*x6 + b1*U2
dtheta/dt     = x4
d^2theta/dt^2 = a2*x2*x6 + b2*U3

What's the problem ?

[Pallov Anand]

I need to write dphid / dt instead of dphi / dt, where

phid(1) = 0;

phid(n+1) = asin(Ux(n)*sin(xd(5,n)) - Uy(n)*cos(xd(5,n)));

Answer 1

[Star Strider]

The diff function to calcualte the derivative is part of the Symbolic Math Toolbox.

To calculate a numerical derivative, use the gradient function.