## 4th order approximations of the second order derivative using different approximation schemes

The aim of this program is to derive a second-order derivative of the fourth-order approximation scheme and then compare its results over a function exp(x)*cos(x). In this program, the first 3 numerical approximations scheme will be derived which are Central difference (CD), Skewed right-sided difference (SRD), and Skewed left-sided difference (SLD) using the Taylor table method. Results will be shown by a comparison plot that will compare the absolute error between the above-mentioned scheme by writing a MatLab program. In last why the skewed scheme is useful is discussed in detail.

1. Central difference (CD) Scheme using Taylor table

The second-order derivative of the fourth-order approximation can be approximated as:-

(del^2f)/(delx^2) approx f''(x) + Deltax^4f^6(x)*Constant = af(i+2) + bf(i+1) + cf(i) + df(i-1) +ef(i-2)                                                                                                          (1)

The Taylor table for the above-mentioned approximation

 Taylor table terms Numerical stencil f(i) ∆xf'(i) ((∆x)2f''(i))/2 ((∆x)3f3(i))/6 ((∆x)4f4(i))/24 ((∆x)5f5(i))/120 ((∆x)6f6(i))/720 af(i+2) a (2)a (2)2 (2)3a (2)4a (2)5a (2)6a bf(i+1) b b b b b b b cf(i) c 0 0 0 0 0 0 df(i-1) d -d d -d d -d d ef(i-2) e (-2)e (-2)2e (-2)3e (-2)4e (-2)5e (-2)6e co-efficient in equation ∑= 0 0 1 0 0 ? Constant

Here coefficient of a second-order derivative term is equal to 1 and the constant term of 6th order derivative is the error respectively.

The equation which need to be solved are:-

a + b + c + d + e = 0   eq(1)

2a + b - d - 2e = 0  eq(2)

4a + b + d + 4e = 1  eq(3)

8a + b - d - 8e = 0  eq(4)

16a + b + d + 16e = 0  eq(5)

32a + b - d - 32e = ?    eq(6)

we have 5 coefficient and first 5 equation will be sufficient to solve it.

2. Skewed right-sided difference (SRD) using Taylor table

The second-order derivative of the fourth-order approximation can be approximated as:-

(del^2f)/(delx^2) approx f''(x) + Deltax^4f^6(x)*Constant = af(i) + bf(i+1) + cf(i+2) + df(i+3) + ef(i+4) + gf(i+5)                                                                                        (2)

The Taylor table for the above mentioned approximation:-

 Taylor table terms Numerical stencil f(i) ∆xf'(i) ((∆x)2f''(i))/2 ((∆x)3f3(i))/6 ((∆x)4f4(i))/24 ((∆x)5f5(i))/120 ((∆x)6f6(i))/720 af(i) a 0 0 0 0 0 0 bf(i+1) b b b b b b b cf(i+2) c 2c 4c 8c 16c 32c 64c df(i+3) d 3d 9d 27d 81d 243d 729d ef(i+4) e 4e 16e 64e 256e 1024e 4096e gf(i+5) g 5g 25g 125g 625g 3125g 15625g co-efficient in equation ∑= 0 0 1 0 0 0 Constant

Here coefficient of a second-order derivative term is equal to 1 and the constant term of 6th order derivative is the error respectively.

The equation which need to be solved are:-

a + b + c + d + e = 0                           eq(1)

b + 2c + 3d +4e + 5g= 0                      eq(2)

b + 4c + 9d + 16e + 25g = 2                eq(3)

b + 8c + 27d + 64e + 125g = 0            eq(4)

b + 16c + 81d +256e + 625g = 0         eq(5)

b + 32c + 243d +1024e + 3125g = 0    eq(6)

3. Skewed left-sided difference (SLD) using Taylor table

The second-order derivative of the fourth-order approximation can be approximated as:-

(del^2f)/(delx^2) approx f''(x) + Deltax^4f^6(x)*Constant = af(i) + bf(i-1) + cf(i-2) + df(i-3) + ef(i-4) + gf(i-5)                                                                                            (3)

The Taylor table for the above mentioned approximation:-

 Taylor table terms Numerical stencil f(i) ∆xf'(i) ((∆x)2f''(i))/2 ((∆x)3f3(i))/6 ((∆x)4f4(i))/24 ((∆x)5f5(i))/120 ((∆x)6f6(i))/720 af(i) a 0 0 0 0 0 0 bf(i-1) b (-1)b b (-1)b b (-1)b b cf(i-2) c (-2)c 4c (-8)c 16c (-32)c 64c df(i-3) d (-3)d 9d (-27)d 81d (-243)d 729d ef(i-4) e (-4)e 16e (-64)e 256e (-1024)e 4096e gf(i-5) e (-5)e 25e (-125)g 625g (-3125)g 15625e co-efficient in equation ∑= 0 0 1 0 0 0 Constant

Here coefficient of a second-order derivative term is equal to 1 and the constant term of 6th order derivative is the error respectively.

The equation which need to be solved are:-

a + b + c + d + e = 0                               eq(1)

-(b + 2c + 3d +4e + 5g)= 0                      eq(2)

b + 4c + 9d + 16e + 25g = 2                    eq(3)

-(b + 8c + 27d + 64e + 125g) = 0            eq(4)

b + 16c + 81d +256e + 625g = 0             eq(5)

-(b + 32c + 243d +1024e + 3125g) = 0    eq(6)

Solving the equations for different scheme:-

1. Central difference (CD) Scheme

(del^2f)/(delx^2) = (-f(i+2) + 16f(i+1) - 30f(i) + 16f(i-1) - f(i-2))/12

2. Skewed right-sided difference (SRD) Scheme

(del^2f)/(delx^2) = (45f(i) - 154f(i+1) + 214f(i+2) - 156f(i+3) + 61f(i+4) - 10f(i+5))/12

3. Skewed left-sided difference (SLD) Scheme

(del^2f)/(delx^2) = (45f(i) - 154f(i+1) + 214f(i+2) - 156f(i+3) + 61f(i+4) - 10f(i+5))/12

After solving for coeffecients of different schemes, it is clearly observed that all the schemes are 4th order accurate as the constant term of Δx^4f^6(x) is not zero and also the absolute value of the overall term with constant is bigger as compared to higher order derivative if calculated.

The Matlab program to get the coefficients is attached below. Here the program solves the equation made for different schemes and finds the coefficients of the respective schemes.

%% Code_Name: coefficients_of_second_order_derivative_for_different_scheme.m
clear all
close all
clc
syms a b c d e g
%% Central Difference
eqn1 = a + b + c + d + e == 0;
eqn2 = 2*a + b - d - 2*e == 0;
eqn3 = 2^2*a + b + d + 2^2*e == 2;
eqn4 = 2^3*a + b - d - 2^3*e == 0;
eqn5 = 2^4*a + b + d + 2^4*e == 0;
[A,B] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5], [a, b, c, d, e]);
X1 = linsolve(A,B);% Solving for coefficient of equation
% Coefficient of Central difference
fprintf('The coefficient of Central Difference')
coeff_a = X1(1)
coeff_b = X1(2)
coeff_c = X1(3)
coeff_d = X1(4)
coeff_e = X1(5)

%% skewed right-sided difference
eqn1 = a + b + c + d + e + g == 0;
eqn2 = b + 2*c+3*d+4*e+5*g == 0;
eqn3 = b + 2^2*c + 3^2*d+ 4^2*e + 5^2*g == 2;
eqn4 = b + 2^3*c + 3^3*d + 4^3*e + 5^3*g == 0;
eqn5 = b + 2^4*c + 3^4*d + 4^4*e + 5^4*g == 0;
eqn6 = b + 2^5*c + 3^5*d + 4^5*e + 5^5*g == 0;
[A,B] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6], [a, b, c, d, e, g]);
X2 = linsolve(A,B);% Solving for coefficient of equation
% Coefficient of Skewed right-sided difference
fprintf('The coefficient of Skewed right-sided Difference')
coeff_a = X2(1)
coeff_b = X2(2)
coeff_c = X2(3)
coeff_d = X2(4)
coeff_e = X2(5)
coeff_g = X2(6)

%% skewed left-sided difference
eqn1 = a + b + c + d + e + g == 0;
eqn2 = -(b + 2*c + 3*d + 4*e + 5*g) == 0;
eqn3 = b + 2^2*c + 3^2*d+ 4^2*e + 5^2*g == 2;
eqn4 = -(b + 2^3*c + 3^3*d + 4^3*e + 5^3*g) == 0;
eqn5 = b + 2^4*c + 3^4*d + 4^4*e + 5^4*g == 0;
eqn6 = -(b + 2^5*c + 3^5*d + 4^5*e + 5^5*g) == 0;
[A,B] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6], [a, b, c, d, e, g]);
X3 = linsolve(A,B);% Solving for coefficient of equation
% Coefficient of Skewed left-sided difference
fprintf('The coefficient of Skewed left-sided Difference')
coeff_a = X3(1)
coeff_b = X3(2)
coeff_c = X3(3)
coeff_d = X3(4)
coeff_e = X3(5)
coeff_g = X3(6)

The a, b and other corresponding values for different schemes are obtained with the help of the above code.

The code is written to evaluate the Second-Order derivative. The code will evaluate the second-order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations schemes derived above. The value has computed over x=pi/4 with different values of dx from pi/40 to pi/40000. The program and its corresponding function are attached below.

%% Code_Name: second_order_derivative_fourth_order_approximation.m
clear all
close all
clc

%% The function that has to be computed is exp(x)cos(x);

% Analytical Function = exp(x)cos(x);
% Analytical Second Order Derivative (ASOD)
% f''(x) = -2*exp(x)*sin(x);

%% Computing value at x = pi/4;
x = pi/4;
% for computing numerically...let us compute by taking a small change in x
dx = linspace(pi/40, pi/40000, 100);% Step Value

%% Computing value both analytically and numerically

% Analytical Second Order Derivative (ASOD)
ASOD = -2*exp(x)*sin(x)

% Numerical Derivative

[CD SRD SLD CD_error SRD_error SLD_error] = sod_foa_fun(x,dx,ASOD)

%% 4. Plotting bar graph

% Showing error plot for Central Difference Scheme
loglog(dx,CD_error,'r');
hold on

% Showing error plot for Skewed right-sided Difference Scheme
loglog(dx,SRD_error,'b');
hold on

% Showing error plot for Skewed left-sided Difference Scheme
loglog(dx,SLD_error,'k');

%% 5. Appropriate legends and label

xlabel('dx')% Step size
ylabel('error')% absolute difference between numerical and analytical value at different step size
legend('CD error','SRD error', 'SLD error')% shows error values legend on loglog scale
title('Error Comparison of Second-Order Derivative')

%% Function_Name: sod_foa_fun.m
function [CD SRD SLD CD_error SRD_error SLD_error] = sod_foa_fun(x,dx,ASOD)

% Analytical Function = exp(x)*cos(x);
% Analytical Second Order Derivative (ASOD)
% f'(x) = -2*exp(x)*sin(x);

% Numerical Derivative
for i=1:length(dx)
%% Central difference (CD) Scheme using Taylor table
%CD = (-f(x+2dx)+16f(x+dx)-30f(x)+16f(x-dx)-f(x-2dx))/(12dx)
CD(i) = (-(exp(x+2*dx(i))*cos(x+2*dx(i)))...
+16*(exp(x+dx(i))*cos(x+dx(i)))...
-30*(exp(x)*cos(x))...
+16*(exp(x-dx(i))*cos(x-dx(i)))...
-(exp(x-2*dx(i))*cos(x-2*dx(i))))/(12*(dx(i))^2)
CD_error(i) = abs(CD(i) - ASOD)

%% Skewed right-sided difference (SRD) using Taylor table
% SRD = (45f(x)-154f(x+dx)+214f(x+2dx)-156f(x+3dx)+61f(x+4dx)-10f(x+5dx)/(12dx)

SRD(i) = (45*(exp(x)*cos(x))-154*(exp(x+dx(i))*cos(x+dx(i)))...
+214*(exp(x+2*dx(i))*cos(x+2*dx(i)))...
-156*(exp(x+3*dx(i))*cos(x+3*dx(i)))...
+61*(exp(x+4*dx(i))*cos(x+4*dx(i)))...
-10*(exp(x+5*dx(i))*cos(x+5*dx(i))))/(12*(dx(i))^2)
SRD_error(i) = abs(SRD(i) - ASOD)

%% Skewed left-sided difference (SLD) using Taylor table
% SRD = (45f(x)-154f(x+dx)+214f(x+2dx)-156f(x+3dx)+61f(x+4dx)-10f(x+5dx)/(12dx)
SLD(i) = (45*(exp(x)*cos(x))...
-154*(exp(x-dx(i))*cos(x-dx(i)))...
+214*(exp(x-2*dx(i))*cos(x-2*dx(i)))...
-156*(exp(x-3*dx(i))*cos(x-3*dx(i)))...
+61*(exp(x-4*dx(i))*cos(x-4*dx(i)))...
-10*(exp(x-5*dx(i))*cos(x-5*dx(i))))/(12*(dx(i))^2)
SLD_error(i) = abs(SLD(i) - ASOD)
end
end

The program will plot a curve that will compare the absolute error between the above-mentioned schemes (for x value of pi/4 and change in dx value from pi/40 to pi/4000) whose image is attached below. It is clearly observed from the plot as the dx value decreases, the respective error of the scheme is also reducing until some value in between e-2 to e-3 and on further decreasing the dx value, the error starts increasing. So a check is required to find so that the results evaluated are correct. Moreover, the error of the Central difference scheme is less as compared to the other two schemes.

The central difference scheme can not be useful in a way as compared to the skewed scheme when we don't have the right-hand side or left-hand side step values of any function. for example, in equation (1) above for central difference scheme if we don't have the values for f(i+2) & f(i+1) or f(i-1) & f(i-2), in that we will not be able to find the result and the central scheme will not be useful at all.