Find the derivative of a piecewise function and plot the result

Hi, I have the following code which includes a piecewise function "s" of theta and I would like to differentiate this function with respect to theta and plot the result after.

 L = 2.5;
 betha1 = 50*pi/180;
 betha2 = 260*pi/180;

 curv_arm = @(theta) L/2*(1-cos((pi*(theta-betha1))/betha1));
 curv_par = @(theta) L*(1-(theta-(100*pi/180))/(betha2)+((1)/(2*pi))*(sin(2*pi*(theta-(100*pi/180))/(betha2))));

 syms theta
 s = @(theta) [0*((0<theta) & (theta<50*pi)) + curv_arm(theta).*((50*pi/180<theta) & (theta<100*pi/180)) + 
 curv_par(theta).*((100*pi/180<theta) & (theta<2*pi))];
 fplot(s, [0,6], 'b')

I have managed to plot the current function but I can't differentiate with the command diff for some reason

Answer 1

[Star Strider]

The Symbolic Math Toolbox diff function will not differentiate it because of the discontinuities, either using my function or your original symbolic function for ‘s’. (The derivative does not exist at the discontinuities.)

The best you can do is to evaluate it numerically over a suitably fine vector for ‘theta’ (use the linspace function) and then use the gradient function to calculate the numerical derivative.

Then you can plot both the function and its numerical derivative as functions of ‘theta’.

[Edgar Dominguez Ugalde]

OK, now my question is: how do I write the arguments in | gradient| once I have defined the | linspace| vector? I am trying with this but I'm getting an error:

v = linspace(0,6);
vel = @(theta) gradient(@(theta) s(theta),v);

I want to plot vel

[Star Strider]

Try this:

L = 2.5;
betha1 = 50*pi/180;
betha2 = 260*pi/180;
curv_arm = @(theta) L/2*(1-cos((pi*(theta-betha1))/betha1));
curv_par = @(theta) L*(1-(theta-(100*pi/180))/(betha2)+((1)/(2*pi)).*(sin(2*pi*(theta-(100*pi/180))/(betha2))))
s = @(theta) [0*((0<theta) & (theta<50*pi)) + curv_arm(theta).*((50*pi/180<theta) & (theta<100*pi/180)) + curv_par(theta).*((100*pi/180<theta) & (theta<2*pi))];

thta = linspace(0, 2*pi, 500);
sv = s(thta);
dsdthta = gradient(sv, mean(diff(thta)));

figure(1)
plot(thta, sv, '-b')
hold on
plot(thta, dsdthta, '--r')
hold off
grid
legend('s(\theta)', 'vel(\theta)')

[Edgar Dominguez Ugalde]

Thanks again, you saved me! :)

[Star Strider]

As always, my pleasure!

Please be careful with the simple ‘less than’ and ‘greater than’. The derivative will appear smoother with:

s = @(theta) [0*((0<=theta) & (theta<50*pi)) + curv_arm(theta).*((50*pi/180<=theta) & (theta<100*pi/180)) + curv_par(theta).*((100*pi/180<=theta) & (theta<=2*pi))];

so the discontinuities are minimised. Experiment to get the result you want.

[Edgar Dominguez Ugalde]

As a last question: now I'm trying to calculate the second and third derivative of the function "s". I looked up for this in the diff documentation and it says that I just need to add the 2 and the 3 after the comma in the command diff for the second and third derivative, respectively but it's not working for me (I'm getting the same plot). What would you suggest for this?

[Star Strider]

You likely cannot take the derivative of a discontinuous function. What you read is correct for the Symbolic Math Toolbox diff function. The core MATLAB diff function would shorten the output vector by 2 and 3 elements respectively.

I would use the numeric gradient function repeatedly to generate the second and third numeric derivatives.

If you want to take the symbolic derivatives, you would have to take the derivatives of the separate segments and then combine them. I do not claim that will be the correct derivative of your entire function, since I have never thought about it before.

The numerical derivatives using the gradient function may also not be appropriate, although it is the only approach I can imagine. That is the problem with discontinuous functions.

[Edgar Dominguez Ugalde]

I'm following your recommendations on using the gradient function to obtain the second and third derivatives, but for some reason the plot I'm getting is the exact same as the first derivative. The code I wrote is:

L = 2.5;
 betha1 = 50*pi/180;
 betha2 = 260*pi/180;
curv_arm = @(theta) L/2*(1-cos((pi*(theta-betha1))/betha1));
curv_par = @(theta) L*(1-(theta-(100*pi/180))/(betha2)+((1)/(2*pi)).*(sin(2*pi*(theta-(100*pi/180))/(betha2))))
s = @(theta) [0*((0<theta) & (theta<50*pi)) + curv_arm(theta).*((50*pi/180<theta) & (theta<100*pi/180)) + curv_par(theta).*((100*pi/180<theta) & (theta<2*pi))];

 %First derivative
thta = linspace(0,2*pi,500);
sv = s(thta);
dsdthta = gradient(sv,mean(diff(thta)));
 %Second derivative
dsdthta2 = gradient(sv,mean(diff(thta,2)));
%Plots
 figure (1)
plot(thta, sv, '-b')
grid
figure (2)
plot(thta, dsdthta,'-r')
grid
figure (3)
plot(thta,dsdthta2,'-g')

And the figure 3 should be the second derivative but it is showing the same plot as figure (2)

[Star Strider]

Not quite. I should have been clearer.

The second derivative should be:

dsdthta2 = gradient(gradient(sv,mean(diff(thta))));

or a bit more efficiently:

dsdthta2 = gradient(dsdthta,mean(diff(thta)));

Use gradient twice to create the second derivative (taking the derivative of the first derivative).

[Edgar Dominguez Ugalde]

Now I get it! Thanks man, you da real MVP

[Star Strider]

Thank you!

(I removed the ‘spam’ warning.)