Scilab: How to use ´numderivative´ function

Question

• I am a new user of Scilab and I am not a mathematician.
• As my end goal, I want to calculate (and plot) the derivative of a piece-wise defined function, see here.
• I tried to start small and just use a simple (continuous) function: f(x) = 3*x.
• My Google-Fu lead me to the numderivative function.
• Problem: It seems that I do not understand how the argument x works since the result is not a 1D-array, instead, it is a matrix.
• Update 1: Maybe I use the wrong function and diff is the way to go. But what is then the purpose of numderivative?

PS: Is this the right place to ask Scilab-related questions? It seems that there are several StackOverflow communities where Scilab-related questions are asked.

---

// Define limits
x0 = 0;
x1 = 2;
// Define array x for which the derivative will be calculated. 
n = 100;
x = linspace (x0, x1, n);

// Define  function f(x)
deff('y=f(x)','y=3*x');

// Calculate derivative of f(x) at the positions x
myDiff = numderivative(f,x)

(I expect the result 3 3 and not a matrix.)

Answer 1

numderivative(f,x) will give you the approximated derivative/Jacobian of f at the single vector x. For your example it yields 3 times the identity matrix, which is the expected result since f(x)=3*x. If you rather need the derivative of f considered as a function of a single scalar variable at x=1 and x=2, then numderivative is not convenient as you would have to make an explicit loop. Just code the formula yourself (here first order formula) :

// Define limits
x0 = 0;
x1 = 2;
// Define array x for which the derivative will be calculated. 
n = 100;
x = linspace (x0, x1, n);

// Define  function f(x)
deff('y=f(x)','y=3*x');x = [1 2];
h = sqrt(%eps);
d = (f(x+h)-f(x))/h;

The formula can be improved (second order or complex step formula).

Answer 2

I accidentally (but fortunately) found the elegant solution to your end goal (and also to mine) of calculating and plotting the derivatives through the SciLab documentation itself via the splin (https://help.scilab.org/docs/6.1.1/en_US/splin.html) and interp (https://help.scilab.org/docs/6.1.1/en_US/interp.html) functions. This elegantly works not only for piece-wise defined functions but also for functions defined numerically by data.

The former will give you the corresponding first derivative of your array values while the latter will give you up to the third derivative. Both are actually meant to be used in conjunction because the interp function requires the corresponding derivatives for your y values, which you can easily get through the splin function. Using this method is better than using the diff function because the output array is of the same size as the original one.

Illustrating this through your code example:

// Define limits
x0 = 0;
x1 = 2;
// Define array x for which the derivatives will be calculated. 
n = 100;
x = linspace (x0, x1, n);

// Define  function f(x)
deff('y=f(x)','y=3*x');

// Calculate derivative of f(x) at the positions x
myDiff = splin(x, f(x));

// Calculate up to the third derivative of f(x) at the positions x
// The first derivatives obtained earlier are required by this function
[yp, yp1, yp2, yp3]=interp(x, x, f(x), myDiff);

// Plot and label all the values in one figure with gridlines
plot(x', [yp', yp1', yp2', yp3'], 'linewidth', 2)
xgrid();

// Put the legend on the upper-left by specifying 2 after the list of legend names
legend(['yp', 'yp1', 'yp2', 'yp3'], 2);