7
$\begingroup$

I wish to compute the derivative of a vector through numerical differentiation. Let's say, we use a standard 2nd order central difference scheme, to arrive at a differentiation matrix, and apply it on a numerical vector, $u$.

$$u' = D u$$

At the top and bottom few elements of the resulting derivative vector, I get spurious numerical round-off even by using complicated stencil.

Now, if I remove the mean from every element of my vector, and then pre-multiply it with my differentiation matrix $D$,

$$\tilde u' = D (u - mean(u))$$

it seems to get better results (i.e. better accuracy), although not fully accurate.

Did I happen to get lucky by chance for this specific vector I am working with, or is there a fundamental mathematical reason why the same differentiation matrix performed remarkably on the mean-removed input vector ?

$\endgroup$
  • $\begingroup$ Are you using the discrete mean or the continuous mean? $\endgroup$ – KyleW Mar 11 '17 at 19:46
  • $\begingroup$ Reading your description I think there's a chance the coefficients in $D$ are somehow wrong. Because it's differentiation, and so rows of $D$ should sum to zero, adding a constant number to every element in a vector should make no real difference. $\endgroup$ – Kirill Mar 11 '17 at 20:05
  • $\begingroup$ @KyleW, the vector to be differentiated is composed of just discrete numerical values. So, I am adding up all their values and dividing by the length of the vector to compute the mean numerically $\endgroup$ – Krishna Mar 11 '17 at 21:35
  • $\begingroup$ @Kirill, yes I realise that. But this was one of the options to my question suggested by one of the members: scicomp.stackexchange.com/a/26358/16900 $\endgroup$ – Krishna Mar 11 '17 at 21:39
  • 2
    $\begingroup$ @Krishna There's nothing spurious about the magnitude of roundoff errors in your question over there. You also really shouldn't test "mean removal" on just constant vectors: I get why it works for the constant vector, but that's not important for finite differencing. Is there a misunderstanding here about how roundoff errors work, and what amount of roundoff is problematic? I fear you might be expecting it to produce exact zeros, or something like that. The numbers you get are consistent with relative perturbations in input data of size $O(\epsilon_{\mathrm{mach}})$. $\endgroup$ – Kirill Mar 12 '17 at 20:49
5
$\begingroup$

I took an example from your other question and tried to demonstrate here what is going on. So, for completeness:

We are applying the 6th order finite-difference scheme for numerical differentiation using the following code. First, let's apply it on the constant function of different scale ($C_1, C_2, C_3$):

$$ U(x)=C,\quad x\in[0,1] $$ Since $\frac{\text{d}U(x)}{\text{d}x}=0$, the expected result of differentiation is a zero-vector.

I used $N$=100 discretization points for

  1. $C_1=2.5545007\cdot10^2$, shown in blue
  2. $C_2=2.5545007\cdot10^4$, shown in red (original scale of the problem you described)
  3. $C_3=2.5545007\cdot10^6$, shown in black

The function (solid lines), as well as it's derivative computed without mean subtraction (circle markers) and with mean subtraction (cross markers) is in the figure below. All crosses are exactly at $0$ (I had to assign a dummy value to display them on the semilog plot), which is an expected behavior for a constant function, as Kirill mentioned in a comment. Because you are passing to the differentiation subroutine (or $U-\text{mean}(U)=0$) vector of zeros. enter image description here

Now, consider a tiny bit more complicated function

$$ U(x)=C\sin(x),\quad x\in[0,1] $$ Since $\frac{\text{d}U(x)}{\text{d}x}=C\cos(x)$, you would also expect to get a zero-vector for $\frac{\text{d}U(x)}{\text{d}x}-C\cos(x)$, and the scales of the problem are similar to the ones above with a constant function. enter image description here

For $\sin$, the mean subtraction did nothing for improving the accuracy of the computed derivative.

So in general, subtracting the mean value (i.e. a constant) does not help too much for functions more complicated than a constant.

Regarding the problem with numerical differentiation being sensitive to round off, you already got several good answers in the comments to both questions.

Matlab code for completeness. The function dss006 is available here.

xmin=0;    xmax=1;
N=100;
x=linspace(xmin,xmax,N);

fudge=1E-17;
cU1=2.5545007E+2;    cU2=2.5545007E+4;    cU3=2.5545007E+6;
U1=cU1*ones(N,1);    U2=cU2*ones(N,1);    U3=cU3*ones(N,1);
W1=cU1*sin(x);    W2=cU2*sin(x);    W3=cU3*sin(x);
dUx1=dss006(xmin,xmax,N,U1);
dUx2=dss006(xmin,xmax,N,U2);
dUx3=dss006(xmin,xmax,N,U3);
dUx1p=dss006(xmin,xmax,N,U1-mean(U1))+fudge;
dUx2p=dss006(xmin,xmax,N,U2-mean(U2))+fudge;
dUx3p=dss006(xmin,xmax,N,U3-mean(U3))+fudge;
dWx1=dss006(xmin,xmax,N,W1);
dWx2=dss006(xmin,xmax,N,W2);
dWx3=dss006(xmin,xmax,N,W3);
dWx1p=dss006(xmin,xmax,N,W1-mean(W1))+fudge;
dWx2p=dss006(xmin,xmax,N,W2-mean(W2))+fudge;
dWx3p=dss006(xmin,xmax,N,W3-mean(W3))+fudge;

%% ONLY PLOTTING AFTER THIS POINT
MarkerSize=7;    LineWid=2.5;    ff=22;  MyLineWidth=1.0;
set(groot,'DefaultTextInterpreter','latex');
set(groot,'defaultLegendInterpreter','latex');
set(groot,'defaultAxesTickLabelInterpreter','latex'); 

figure(1)
ph(1)=semilogy(x,abs(U1),'b-','LineWidth',MyLineWidth);hold on;
p1_str='$|U_1(x)|,\;C_1=2.5545007\cdot 10^2$';
ph(2)=semilogy(x,abs(U2),'r-','LineWidth',MyLineWidth);hold on;
p2_str='$|U_2(x)|,\;C_2=2.5545007\cdot 10^4$';
ph(3)=semilogy(x,abs(U3),'k-','LineWidth',MyLineWidth);hold on;
p3_str='$|U_3(x)|,\;C_3=2.5545007\cdot 10^6$';
ph(4)=semilogy(x,abs(dUx1),'bo','LineWidth',MyLineWidth);hold on;
p4_str='$|\widetilde{dU_1/dx}(x)-0|$';
ph(5)=semilogy(x,abs(dUx2),'ro','LineWidth',MyLineWidth);hold on;
p5_str='$|\widetilde{dU_2/dx}(x)-0|$';
ph(6)=semilogy(x,abs(dUx3),'ko','LineWidth',MyLineWidth);hold on;
p6_str='$|\widetilde{dU_3/dx}(x)-0|$';
ph(7)=semilogy(x,abs(dUx1p),'bx','LineWidth',MyLineWidth);hold on;
p7_str='$|\widetilde{dU_1/dx}(x)-0|$,$-$mean($U_1$)';
ph(8)=semilogy(x,abs(dUx2p),'rx','LineWidth',MyLineWidth);hold on;
p8_str='$|\widetilde{dU_2/dx}(x)-0|$,$-$mean($U_2$)';
ph(9)=semilogy(x,abs(dUx3p),'kx','LineWidth',MyLineWidth);hold on;
p9_str='$|\widetilde{dU_3/dx}(x)-0|$,$-$mean($U_3$)';
title('$U(x)=C$');    grid on; set(gca,'FontSize',ff-2);
xlabel('$x$','FontSize',ff); 
ylabel('$|U(x)|, |\widetilde{\frac{dU}{dx}}(x)-\frac{dU}{dx}(x)|$','FontSize',ff);
h=legend([ph(1) ph(2) ph(3) ph(4) ph(5) ph(6) ph(7) ph(8) ph(9)],p1_str,p2_str,p3_str,p4_str,p5_str,p6_str,p7_str,p8_str,p9_str,'Location','Best');
h.FontSize=ff-6;    h.EdgeColor=[1. 1. 1.];
ylim([1E-17 1E+9]);
set(gca,'YTick',[1E-15 1E-12 1E-9 1E-6 1E-3 1E-0 1E+3 1E+6 1E+9]);

figure(2)
ph(1)=semilogy(x,abs(W1),'b-','LineWidth',MyLineWidth);hold on;
p1_str='$|U_1(x)=C_1\sin(x)|,\;C_1=2.5545007\cdot 10^2$';
ph(2)=semilogy(x,abs(W2),'r-','LineWidth',MyLineWidth);hold on;
p2_str='$|U_2(x)=C_2\sin(x)|,\;C_2=2.5545007\cdot 10^4$';
ph(3)=semilogy(x,abs(W3),'k-','LineWidth',MyLineWidth);hold on;
p3_str='$|U_3(x)=C_3\sin(x)|,\;C_3=2.5545007\cdot 10^6$';
ph(4)=semilogy(x,abs(dWx1-cU1*cos(x)),'bo','LineWidth',MyLineWidth);hold on;
p4_str='$|\widetilde{dU_1/dx}(x)-C_1\cos(x)|$';
ph(5)=semilogy(x,abs(dWx2-cU2*cos(x)),'ro','LineWidth',MyLineWidth);hold on;
p5_str='$|\widetilde{dU_2/dx}(x)-C_2\cos(x)|$';
ph(6)=semilogy(x,abs(dWx3-cU3*cos(x)),'ko','LineWidth',MyLineWidth);hold on;
p6_str='$|\widetilde{dU_3/dx}(x)-C_3\cos(x)|$';
ph(7)=semilogy(x,abs(dWx1p-cU1*cos(x)),'bx','LineWidth',MyLineWidth);hold on;
p7_str='$|\widetilde{dU_1/dx}(x)-C_1\cos(x)|$,$-$mean($U_1$)';
ph(8)=semilogy(x,abs(dWx2p-cU2*cos(x)),'rx','LineWidth',MyLineWidth);hold on;
p8_str='$|\widetilde{dU_2/dx}(x)-C_2\cos(x)|$,$-$mean($U_2$)';
ph(9)=semilogy(x,abs(dWx3p-cU3*cos(x)),'kx','LineWidth',MyLineWidth);hold on;
p9_str='$|\widetilde{dU_3/dx}(x)-C_3\cos(x)|$,$-$mean($U_3$)';
title('$U(x)=C\sin(x)$');    grid on; set(gca,'FontSize',ff-2);
xlabel('$x$','FontSize',ff);
ylabel('$|U(x)|, |\widetilde{\frac{dU}{dx}}(x)-\frac{dU}{dx}(x)|$','FontSize',ff);    
h=legend([ph(1) ph(2) ph(3) ph(4) ph(5) ph(6) ph(7) ph(8) ph(9)],p1_str,p2_str,p3_str,p4_str,p5_str,p6_str,p7_str,p8_str,p9_str,'Location','Best');
h.FontSize=ff-6;    h.EdgeColor=[1. 1. 1.];
ylim([1E-17 1E+9]);
set(gca,'YTick',[1E-15 1E-12 1E-9 1E-6 1E-3 1E-0 1E+3 1E+6 1E+9]);
$\endgroup$
  • $\begingroup$ Much appreciate your complete answer $\endgroup$ – Krishna Apr 15 '17 at 20:23

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.