Calculating a derivative in a loop?

Question

I am not sure where else to post this, so here goes.

I have a velocity profile ($v_{\phi}$) for a radial distance ($r$). From this velocity profile I am trying calculate the shear rate ($y$) which can be calculated by ( Eq 3 in Image): $$ y = \frac{\partial v_{\phi}}{ \partial r} - \frac{v_{\phi}}{r} = r \frac{\partial}{\partial r} \left( \frac{v_{\phi}}{ r}\right ) $$ Please find the image of the equation from the papers attached,

Screenshot of the equation for the paper for your kind consideration.

I have implemented the equation as such,

for i = 1: size(v_phi,1)

    y (i) = ( r(i) ) * (( ( v_phi (i) / r(i) ) - ( v_phi (i+1) / r(i+1) ) ) / ( r(i+1) - r(i)) );

    y1 (i) = ( ( ( v_phi(i) / r(i) ) - ( v_phi(i+1) / r(i+1) ) ) * ( v_phi (i) / r(i) ) );


end

However the y and y1 values calculated are very different. I am just wondering, whether I am doing something wrong in the implementation and looking for help if possible.

Many thanks.

Data >>>>

v_phi = [36.99 37.00 35.78 31.43 26.91 23.32 19.46 16.97 14.49 12.12 10.42 8.52 7.25 6.06 5.11 3.92 3.30 2.54 2.11 1.42 1.46 1.46 0.97 0.63 0.63 0.51 0.40 0.91 0.24 0.36 0.63 0.71 0.45 0.26 0.74 ];

r = [3.59 3.76 3.94 4.12 4.29 4.47 4.65 4.82 5.00 5.18 5.35 5.53 5.71 5.88 6.06 6.24 6.41 6.59 6.76 6.94 7.12 7.29 7.47 7.65 7.82 8.00 8.18 8.35 8.53 8.71 8.88 9.06 9.24 9.41 9.50 ];

Answer 1

Your implementation is:

for i = 1: size(v_phi,1)
  y(i) = r(i)*((v_phi(i)/r(i)  - v_phi(i+1)/r(i+1)) / (r(i+1) - r(i)));
  y1(i) = (v_phi(i)/r(i)  - v_phi(i+1)/r(i+1)) * (v_phi(i)/r(i));
end

That means that the discretization you are using is $$ y_i = r_i \left(\frac{v_{\phi,i}}{r_i} - \frac{v_{\phi, i+1}}{r_{i+1}}\right) \left(\frac{1}{r_{i+1} - r_i}\right) \quad \equiv \quad y = r \frac{d}{d r}\left(\frac{v_\phi}{r}\right) $$ and $$ y_{1,i} = \left(\frac{v_{\phi,i}}{r_i} - \frac{v_{\phi, i+1}}{r_{i+1}}\right) \left(\frac{v_{\phi,i}}{r_i}\right) \quad \equiv \quad y_1 = \frac{d v_\phi}{d r} - \frac{v_\phi}{r} $$ From the above we notice that the discretization may be inaccurate and there is at least one bug in your implementation.

Also, the ratio of the two values is $$ \frac{y_i}{y_{1,i}} = \frac{r_i^2}{v_{\phi,i} (r_{i+1} - r_i)} $$ There is no reason why this ratio should be 1, indicating that there is an error in your implementation.

To fix the finite difference accuracy, you have to go to the definition of the discretization scheme. In your case you appear to be using a forward difference scheme. For that situation, as $h \rightarrow 0$, $$ \frac{df(x)}{dx} = \frac{f(x+h) - f(x)}{h} + O(h) $$ The accuracy of your solution depends on the value of $h$. To increase accuracy you should use a central difference scheme at the minimum.

Answer 2

Biswajit Banerjee already indicated that your code has bugs, which is why your results are so far off. I will explain in more detail where these bugs are:

v_phi is a row vector, so size(v_phi, 1) always returns 1 because there is only one row. You want size(v_phi, 2) or better yet, numel(v_phi).

Because you're using forward differences to compute the derivatives, you can only iterate up to the second-to-last element in each of the input vectors, so the loop declaration should start like:

for i = 1:size(v_phi, 2) - 1

After that, you mix up the order of subtraction at least once or twice. That is, where you should be doing y(i+1) - y(i), you instead do y(i) - y(i+1), for example.

Also, as a matter of style, loops can often be replaced by vector operations in Matlab. For example, the following one-liner computes $\frac{dv_{\phi}}{dr}$ over the whole data set:

dvdr = diff(v_phi) ./ diff(r);

Your equation can actually be solved without using any loops.