MATLAB: First and Second Order Central Difference

MATLAB

The 1st order central difference (OCD) algorithm approximates the first derivative according to ,
and the 2nd order OCD algorithm approximates the second derivative according to .
In both of these formulae is the distance between neighbouring x values on the discretized domain.
a.
Write a script which takes the values of the function for and make use of the 1st and 2nd order algorithms to numerically find the values of and . You may use the analytical value of to find initial condtions if required.
b.
Plot your results on two graphs over the range , comparing the analytical and numerical values for each of the derivatives.
c.
Compare each numerical algorithms results by finding the largest value of the relative error between the analytical and numerical results.
Can someone please help with this question? I'm stuck on where to begin really. Thanks! This is what I have so far, but comes up with errors.
clear all
f=@(x) cosh(x)
x=linspace(-4,4,9)
n=length(x)
i=1:n
h=x(i)-x(i-1)
xCentral=x(2:end-1);
dFCentral=(F(i+1)-F(i))/(h);

Best Answer

For starters, the formula given for the first derivative is the FORWARD difference formula, not a CENTRAL difference.
(here, dt = h)
Second: you cannot calculate the central difference for element i, or element n, since central difference formula references element both i+1 and i-1, so your range of i needs to be from i=2:n-1.
f = @(x) cosh(x);
h = 1;
x = -4:h:4;   % this way, you define the desired step size, h, and use it to calculate the x vector
              % to change the resolution, simply change the value of h
% x = linspace(-4,4,9);
n = length(x);
y=f(x);
dy = zeros(n,1); % preallocate derivative vectors
ddy = zeros(n,1);
for i=2:n-1
    dy(i) = (y(i-1)+y(i+1))/2/h;
    ddy(i) = (y(i+1)-2*y(i)+y(i-1))/h^2;
end
% Now when you plot the derivatives, skip the first and the last point
figure;
plot(x,y,'r');
hold on;
plot(x(2:end-1), dy(2:end-1),'b');
plot(x(2:end-1), ddy(2:end-1), 'g');
grid on;
legend('y', 'dy', 'ddy')
Try making h smaller to see how it effects the result. (h=0.1, h=0.01)
Another change you might consider, in order to fill in the first and last point in the derivative is:
for i=1:n
    switch i
        case 1
            % use FORWARD difference here for the first point
            dy(i) = ...
            ddy(i) = ...
        case n
            % use BACKWARD difference here for the last point
            dy(i) = ...
            ddy(i) = ...
        otherwise
            % use CENTRAL difference
            dy(i) = ...
            ddy(i) = ...
    end
end
% Now you can plot all points in the vectors (from 1:n)
Related Question