Strange graph artifact when performing diff() of a high order

I think it has something to do with values going below eps within the code but I'm not sure. Is there a simple fix or do I have to rework the code entirely?

clearvars;

x1 = 0; %Start
x2 = 4*pi; %End
I = 1000; %Steps
h = (x2-x1)/I; %Step size
N = 5; %Number of derivatives

x = linspace(x1,x2,I);
y = exp(1/2*x); %Base function
dydx = zeros(N, I-N); %Preallocation

for n = 1:N
    diffVector = diff(y,n)/h^n; %Calculates the nth derivative
    dydx(n,:) = diffVector(1:I-N); %Stores it in the derivative array
end

plot(x(:,1:I-N),dydx,x,y) %Plots all derivatives

The code is just an exercise to find the nth derivative of a function (in this case exp(1/2x))

When N = 5 it seems to work fine. But if I change it to N = 7 I get some freaky graph output. This may also occur if I change I by a certain amount. Here are the graphs of each:

<<

>>

Answer 1

[dpb]

diff is the simple differential between adjacent points; when you begin subtracting values one after another, loss of precision begins to happen and become noticeable as you've uncovered.

Increasing I exacerbates the problem as the spacing between points decreases with increasing length of the vector over a fixed range because it means the difference between one value and the next in the original vector is smaller so the difference starts out smaller, sorta' like you'd already taken one difference already.

There's no fix numerically easily; only increased precision to not lose bits will let you simply diff results indefinitely. It would be the place to use analytic forms or symbolic or other approximations rather than brute force.

Answer 2

[Roger Stafford]

For the higher derivatives you need to increase the size of h which means decrease the number of points, I. There has to be a balance here between the error due to inaccuracies of the higher order differences on the one hand and the variation of the function throughout the range of n*h on the other hand. The higher the derivative, the larger the range of n+1 points has to be.