# First and Second Order Central Difference

2,438 views (last 30 days)
Eleanor Lee on 3 Dec 2019
Edited: Sushant Powar on 14 Oct 2020
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);

Jim Riggs on 3 Dec 2019
Edited: Jim Riggs on 3 Dec 2019
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)
##### 2 CommentsShowHide 1 older comment
Sushant Powar on 14 Oct 2020
Hi Jim, Just heads up you have got the wrong sign in the following line of code:
dy(i) = (y(i-1)+y(i+1))/2/h;
It should be
dy(i) = (y(i+1)-y(i-1))/2/h;
Cheers!