Plotting differentiation error in matlab

Question

I wrote a matlab code that approximates the derivative of cos(x) using the central difference.

meaning cos(x)'=(cos(x+h)-cos(x-h))/2h approximately.

as h goes from 0.1,0.01,0.001 and so on until 10^-8

This is my code

function [pos,neg,D]=shifted_cos(x)
    for i=1:8
        var=rand(1)*0.5*10^(-5);
        pos(i)=cos(x+10^(-i))+var;
    end
    for j=1:8
        var=rand(1)*0.5*10^(-5);
        neg(j)=cos(x-10^(-j))+var;
    end
    D=pos-neg;
    for k=1:8
        D(k)=(D(k)/(2*10^(-k)));
    end
end

Please note that the variable named var is just random noise that i added, since any "real" system has noise and its not 100% reliable, but i was given that the noise is no larger than 0.5*10^-5

Here is my problem: I am now trying to plot the error as a function of h.

By error i mean the real value of the derivative, minus the value in D (please note that D is a vector)

So I wrote the following lines in the console:

[p,n,d]=shifted_cos(1.2)
h=[0.1,0.01,0.001,0.0001,0.00001,0.000001,0.0000001,0.00000001]
plot(h,abs(-sin(1.2)-d))

and the output im getting is very weird. im seeing a graph, but its just a vertical line.

the values of d are

d =

  Columns 1 through 3

  -0.930475812604331  -0.932134403564042  -0.932748001778061

  Columns 4 through 6

  -0.944125866359780  -0.991297975178052   0.416450071288876

  Columns 7 through 8

   8.360791447226124  10.313974302400553

So it shouldnt be a vertical line...

Can someone shed some light on my problem?

Note: I was asked to plot the error as a function of h at the point x=1.2, and also note that the derivative of cos(x) is -sin(x)

mbschenkel mbschenkel · Accepted Answer · 2014-05-27T19:01:16

I've run just the following code and that does not give me a vertical line. However since h is changing in orders of magnitude, a (double) logarithmic plot will give you much more insight.

d = [  -0.930475812604331  -0.932134403564042  -0.932748001778061 ...
       -0.944125866359780  -0.991297975178052   0.416450071288876 ...
        8.360791447226124  10.313974302400553 ];
h = 10.^-(1:8);

%// plot 1: your plot but with adjusted axes
subplot(1, 3, 1)
plot(h, abs(-sin(1.2)-d), '.-')
xlim([-.01, .2])
ylim([-1, 18])

%// plot 2: your data, just linearly plotted
subplot(1, 3, 2)
plot(abs(-sin(1.2)-d), '.-')

%// plot 3: double logarithmic plot
subplot(1, 3, 3)
loglog(abs(-sin(1.2)-d), '.-')

Plot

PS and not related to your question as such: Have a look here on how to avoid for loops in Matlab and make your code both faster and easier to read (And write).

Plotting differentiation error in matlab

1 Answers