Finding Numerical Derivative using a shift matrix in MATLAB

Question

We were given an assignment to find the numerical first and second derivative using given steps:

1) Define two arrays, x and y=f (x) (Use any function)

2) Define the differential operator as a matrix (e.g. for f' , define the matrix D= (1/2h)(U-L) where U is an Uppershift Matrix, and L is Lowershift Matrix.

3) Calculate the derivative by multiplying the differential operator matrix by the function array.

I have written this code:

%QUESTION 3- DIFFERENTIAL OPERATOR
h=2;
x = 2:h:8
y = x.^2 %the chosen function


n=length(x);
shift = 1; 
U = diag(ones(n-abs(shift),1),shift);
s = -1;
L= diag(ones(n-abs(shift),1),s);
% % the code above creates the upper and lower shift matrix
%
D= ((U-L))/(2*h) %differential operator 
d = y*D %approximates the first derivative of each element in vector -x
d2 = d*D %approximates the second derivative of each element in vector -x

for creating the shift matrix i used this post

Now I got this as a solution:

>> mat2q3

x =

     2     4     6     8


y =

     4    16    36    64


D =

         0    0.2500         0         0
   -0.2500         0    0.2500         0
         0   -0.2500         0    0.2500
         0         0   -0.2500         0


d =

    -4    -8   -12     9


d2 =

    2.0000    2.0000   -4.2500   -3.0000

Please tell me where I went wrong, or how can I improve it.

Answer 1

This code for the most part does what it intends to do: compute the numerical approximation to the derivative. However, the way you're computing the derivative is slightly incorrect. Basically, for each point in your array, you want to subtract the point to the left with the point to the right and divide by 2*h. If you want to do that, you need to post-multiply the vector y and transform it into a column vector.

Transpose y then take D and multiply this with y:

>> d = D*y.'

d =

     4
     8
    12
    -9

However, I'd like to point out that the first and last entry don't make any sense because of your D matrix:

D =

         0    0.2500         0         0
   -0.2500         0    0.2500         0
         0   -0.2500         0    0.2500
         0         0   -0.2500         0

What you are computing here is the central difference for the numerical derivative, which requires a point to the left and a point to the right of the point you are evaluating in your array. Basically, for the first point, you only have the point to the right where the left point is undefined as you are going out of bounds. Similarly for the last point, you only have a point to the left and the point to the right is also out of bounds as it doesn't exist.

The numerical derivative matrix here assumes that the points going out of bounds are 0. It just turns out that you get the derivative correct at this point which is a fluke (i.e. 2*x at x = 1 is 2). The reason why is because the point at x = 0 does give the output of x.^2 being equal to 0, so the omission of this point in the array acts as if it was in your point array to begin with.

The same thing is required for the second derivative:

>> d2 = D*d

d2 =

    2.0000
    2.0000
   -4.2500
   -3.0000

However, bear in mind that for your previous first derivative result d, the last entry is garbage and so if you computed the derivative of this result, you would be using the incorrect right most term when you look at the third entry of d, and so the last two and also the first two entries of d2 are garbage. It just so happens that the second entry of d2 is correct. This is attributed to x = 0 not actually existing in the previous result, but with the assumption that values out of bounds in your point array are 0, we do get the right result.

You should try this on a longer sequence. For example, try doing it on this:

h = 2;
x = 2:h:20;
y = x.^2;

This defines a sequence going from 2 to 20 in steps of h = 2. This is what we get for x, and d:

x =

    2     4     6     8    10    12    14    16    18    20

d =

     4
     8
    12
    16
    20
    24
    28
    32
    36
   -81

As you can see the first and last entry are meaningless as we don't have a left most point for the first entry and a right most point for the second entry to consider. However, the first point by fluke gives us the correct result as we talked about earlier. In general, the first and last entry should not give you the correct results as we don't have a left or right point to account for the derivative. For the rest of the points, you can see that we are correctly computing the derivative, or 2*x.

Let's take a look at the second derivative, d2:

d2 =

    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
  -28.2500
   -9.0000

It so happens that the first two elements are correct, even though they theoretically shouldn't, but we got it by fluke. However, the last two elements aren't correct.

---

What you need to take away from this is to remember one thing. Remember that when taking numerical derivatives with the central difference, you have to trim off the first n elements and the last n elements, where n is the order of the derivative you are considering.