Taylor Table for fourth order accurate finite difference approximation using Octave/MATLAB

Objective:- To compare truncation error values between Central difference scheme, Skewed right-sided difference scheme & Skewed left-sided difference scheme for 4th order approximations of the second-order derivative.

Theory:-  A Taylor series is a series expansion of a function about a point. Taylors series for a function `f(x)` at point about a point `x=a` can be given as-

`f(x)=f(a)+f'(a)(x-a)+(f''(a)(x-a)^2)/(2!)+(f'''(a)(x-a)^3)/(3!)+......+(f^n'(a)(x-a)^n)/(n!)....`

1. Central differencing scheme:- for this type of scheme number of nodes can be calculated by following formula.

Number of nodes = p + q - 1

Where, p = order of derivative = 2

           q = order of approximation = 4.

`therefore` Number of nodes = 2+4-1 = 5

`Delta x^2 (d^2u)/dx^2=af(i-2)+bf'(i-1)+cf''(i)+df'''(i+1)+ef^(iv)(i+2)`

Lets prepare Taylors table for above scheme.

Hence we get 5 equations from above table,

`a+b+c+ d+e=0`

`-2a-b+0c+d+2e=0`

`2a-0.5b+0c+0.5d+2e=1`

`-(4/3)a-(1/6)b+0c+(1/6)d+(4/3)e=0`

`(2/3)a+(1/24)b+0c+(1/24)d+(2/3)e=0`

Lets solve above equation using the matrix inversion method and calculate values for coefficients 'a', 'b', 'c', 'd' & 'e' by creating MATLAB code.

clear all
close all
clc

%Objective is to find unknown coefficient for central differencing scheme
%Creating matrix A of coefficients
A = [1 1 1 1 1; -2 -1 0 1 2; 2 0.5 0 0.5 2; -4/3 -1/6 0 1/6 4/3; 2/3 1/24 0 1/24 2/3]

%Matrix B
B = [0; 0; 1; 0; 0]

%Matrix X consist unknown coefficients 
%Matrix X = [a; b; c; d; e]

%Calculate Matrix X by X = AB
X = AB

Hence, we get a=-0.083333; b = 1.3333; c = -2.5; d = 1.3333; e = -0.083333

`therefore` a = e & b = d

`Delta x^2 (d^2u)/dx^2=af(i-2)-bf'(i-1)+cf''(i)+df^iv(i+1)+ef^v(i+2)`

`therefore` `Delta x^2 (d^2u)/dx^2 =` `(a+b+c+d+e).f(i)+``(-2a-b+0c+d+2e)(f'(i) Deltax)+``(2a+0.5b+0c+0.5d+2e)(f''(i) Deltax^2)+``(-(4/3)a-(1/6)b+0c+(1/6)d+(4/3)e)(f'''(i) Deltax^3)+``((2/3)a+(1/24)b+0c+(1/24)d+(2/3)e)(f'^(iv)(i) Deltax^4)+``(-(32/120)a-(1/120)b+0c+(1/120)d+(32/120)e)(f'^(v)(i) Deltax^5)+``((64/720)a+(1/720)b+0c+(1/720)d+(64/720)e)(f'^(vi)(i) Deltax^6)+.....`

For above equation as 

`a+b+c+ d+e=0`

`-2a-b+0c+d+2e=0`

`-(4/3)a-(1/6)b+0c+(1/6)d+(4/3)e=0`

`(2/3)a+(1/24)b+0c+(1/24)d+(2/3)e=0`

and by putting values of 'a', 'b', 'c', 'd' & 'e' we get `(-(32/120)a-(1/120)b+0c+(1/120)d+(32/120)e)=0`

`therefore` `Delta x^2 (d^2u)/dx^2 =` `(2a+0.5b+0c+0.5d+2e)(f''(i) Deltax^2)+``((64/720)a+(1/720)b+0c+(1/720)d+(64/720)e)(f'^(vi)(i) Deltax^6)+.....`

Now divide above equation by `Delta x^2` on both side

Hence we get,

`(d^2u)/dx^2 =` `(2a+0.5b+0c+0.5d+2e)f''(i)+``((64/720)a+(1/720)b+0c+(1/720)d+(64/720)e)(f'^(vi)(i) Deltax^4)+.....`

Which contains `Deltax^4` i.e. 4th order and higher (above 4) order terms.

2. Skewed right sided differencing scheme:- for this type of scheme number of nodes can be calculated by following formula.

Number of nodes = p + q

Where, p = order of derivative = 2

           q = order of approximation = 4.

`therefore` Number of nodes = 2+4 = 6

`Delta x^2 (d^2u)/dx^2=af(i)+bf'(i+1)+cf''(i+2)+df'''(i+3)+ef^(iv)(i+4)+gf^v(i+5)`

Lets prepare Taylors table for above scheme.

Hence we get 6 equations from above table,

`a+b+c+ d+e+g=0`

`0a+b+2c+3d+4e+5g=0`

`0a+0.5b+2c+(9/2)d+8e+(25/2)g=1`

`0a+(1/6)b+(4/3)c+(27/6)d+(64/6)e+(125/6)g=0`

`0a+(1/24)b+(2/3)c+(81/24)d+(256/24)e+(625/24)g=0`

`0a+(1/120)b+(32/120)c+(243/120)d+(1024/120)e+(3125/120)g=0`

Lets solve above equation using the matrix inversion method and calculate values for coefficients 'a', 'b', 'c', 'd', 'e' & 'g' by creating MATLAB code.

clear all
close all
clc

%Objective is to find unknown coefficient for skewed right side differencing scheme
%Creating matrix A of coefficients
A = [1 1 1 1 1 1; 0 1 2 3 4 5; 0 0.5 2 9/2 8 25/2; 0 1/6 4/3 27/6 64/6 125/6; 0 1/24 2/3 81/24 256/24 625/24; 0 1/120 32/120 243/120 1024/120 3125/120]

%Matrix B
B = [0; 0; 1; 0; 0; 0]

%Matrix Y consist unknown coefficients 
%Matrix Y = [a; b; c; d; e; g]

%Calculate Matrix Y by Y = AB
Y = AB

%for calulating coefficient of(Deltax^6). f^(vi)(i)
z=(Y(1)*0)+(Y(2)*(1/720))+(Y(3)*(64/720))+(Y(4)*(729/720))+(Y(5)*(4096/720))+(Y(6)*(15625/720))

Hence, we get a = 3.75; b = -12.833; c = 17.833; d = -13; e = 5.0833; g = -0.83333

`Delta x^2 (d^2u)/dx^2=af(i)-bf'(i+1)+cf''(i+2)+df'''(i+3)+ef^(iv)(i+4)+gf^v(i+5)`

`therefore` `Delta x^2 (d^2u)/dx^2 =` `(a+b+c+ d+e+g).f(i)+``(0a+b+2c+3d+4e+5g)(f'(i) Deltax)+``(0a+0.5b+2c+(9/2)d+8e+(25/2)g)(f''(i) Deltax^2)+``(0a+(1/6)b+(4/3)c+(27/6)d+(64/6)e+(125/6)g)(f'''(i) Deltax^3)+``(0a+(1/24)b+(2/3)c+(81/24)d+(256/24)e+(625/24)g)(f'^(iv)(i) Deltax^4)+``(0a+(1/120)b+(32/120)c+(243/120)d+(1024/120)e+(3125/120)g)(f'^(v)(i) Deltax^5)+``(0a+(1/720)b+(64/720)c+(729/720)d+(4096/720)e+(15625/720)g)(f^(vi)(i)Deltax^6)+....`

For above equation as  

 `a+b+c+ d+e+g=0`

`0a+b+2c+3d+4e+5g=0`

`0a+0.5b+2c+(9/2)d+8e+(25/2)g=1`

`0a+(1/6)b+(4/3)c+(27/6)d+(64/6)e+(125/6)g=0`

`0a+(1/24)b+(2/3)c+(81/24)d+(256/24)e+(625/24)g=0`

`0a+(1/120)b+(32/120)c+(243/120)d+(1024/120)e+(3125/120)g=0`

and by putting values of 'a', 'b', 'c', 'd', 'e' & 'g' we get `(0a+(1/720)b+(64/720)c+(729/720)d+(4096/720)e+(15625/720)g) = -0.7611` which is not equal to zero.

`therefore` `Delta x^2 (d^2u)/dx^2 =` `(0a+0.5b+2c+(9/2)d+8e+(25/2)g)(f''(i) Deltax^2)+``(0a+(1/720)b+(64/720)c+(729/720)d+(4096/720)e+(15625/720)g)(f'^(vi)(i) Deltax^6)+.....`

Now divide above equation by `Delta x^2` on both side

Hence we get,

`(d^2u)/dx^2 =` `(0a+0.5b+2c+(9/2)d+8e+(25/2)g)f''(i)+``(0a+(1/720)b+(64/720)c+(729/720)d+(4096/720)e+(15625/720)g)(f'^(vi)(i) Deltax^4)+.....`

Which contains `Deltax^4` i.e. 4th order and higher (above 4) order terms.

3. Skewed left sided differencing scheme:- for this type of scheme number of nodes can be calculated by following formula.

Number of nodes = p + q

Where, p = order of derivative = 2

           q = order of approximation = 4.

`therefore` Number of nodes = 2+4 = 6

`Delta x^2 (d^2u)/dx^2=af(i)+bf'(i-1)+cf''(i-2)+df'''(i-3)+ef^(iv)(i-4)+gf^v(i-5)`

Lets prepare Taylors table for above scheme.

Hence we get 6 equations from above table,

`a+b+c+ d+e+g=0`

`0a-b-2c-3d-4e-5g=0`

`0a+0.5b+2c+(9/2)d+8e+(25/2)g=1`

`0a-(1/6)b-(4/3)c-(27/6)d-(64/6)e-(125/6)g=0`

`0a+(1/24)b+(2/3)c+(81/24)d+(256/24)e+(625/24)g=0`

`0a-(1/120)b-(32/120)c-(243/120)d-(1024/120)e-(3125/120)g=0`

Lets solve above equation using the matrix inversion method and calculate values for coefficients 'a', 'b', 'c', 'd', 'e' & 'g' by creating MATLAB code.

clear all
close all
clc

%Objective is to find unknown coefficient for skewed left side differencing scheme
%Creating matrix A of coefficients
A = [1 1 1 1 1 1; 0 -1 -2 -3 -4 -5; 0 0.5 2 9/2 8 25/2; 0 -1/6 -4/3 -27/6 -64/6 -125/6; 0 1/24 2/3 81/24 256/24 625/24; 0 -1/120 -32/120 -243/120 -1024/120 -3125/120]

%Matrix B
B = [0; 0; 1; 0; 0; 0]

%Matrix Z consist unknown coefficients 
%Matrix Z = [a; b; c; d; e; g]

%Calculate Matrix Z by Z = AB
Z = AB

%for calulating coefficient of(Deltax^6). f^(vi)(i)
z=(Z(1)*0)+(Z(2)*(1/720))+(Z(3)*(64/720))+(Z(4)*(729/720))+(Z(5)*(4096/720))+(Z(6)*(15625/720))

Hence, we get a = 3.75; b = -12.833; c = 17.833; d = -13; e = 5.0833; g = -0.83333

`Delta x^2 (d^2u)/dx^2=af(i)-bf'(i+1)+cf''(i+2)+df'''(i+3)+ef^(iv)(i+4)+gf^v(i+5)`

`therefore` `Delta x^2 (d^2u)/dx^2 =` `(a+b+c+ d+e+g).f(i)-``(0a+b+2c+3d+4e+5g)(f'(i) Deltax)+``(0a+0.5b+2c+(9/2)d+8e+(25/2)g)(f''(i) Deltax^2)-``(0a+(1/6)b+(4/3)c+(27/6)d+(64/6)e+(125/6)g)(f'''(i) Deltax^3)+``(0a+(1/24)b+(2/3)c+(81/24)d+(256/24)e+(625/24)g)(f'^(iv)(i) Deltax^4)-``(0a+(1/120)b+(32/120)c+(243/120)d+(1024/120)e+(3125/120)g)(f'^(v)(i) Deltax^5)+``(0a+(1/720)b+(64/720)c+(729/720)d+(4096/720)e+(15625/720)g)(f^(vi)(i)Deltax^6)-....`

For above equation as 

`a+b+c+ d+e+g=0`

`0a-b-2c-3d-4e-5g=0`

`0a+0.5b+2c+(9/2)d+8e+(25/2)g=1`

`0a-(1/6)b-(4/3)c-(27/6)d-(64/6)e-(125/6)g=0`

`0a+(1/24)b+(2/3)c+(81/24)d+(256/24)e+(625/24)g=0`

`0a-(1/120)b-(32/120)c-(243/120)d-(1024/120)e-(3125/120)g=0`

and by putting values of 'a', 'b', 'c', 'd', 'e' & 'g' we get `(0a+(1/720)b+(64/720)c+(729/720)d+(4096/720)e+(15625/720)g) = -0.7611` which is not equal to zero.

`therefore` `Delta x^2 (d^2u)/dx^2 =` `(0a+0.5b+2c+(9/2)d+8e+(25/2)g)(f''(i) Deltax^2)+``(0a+(1/720)b+(64/720)c+(729/720)d+(4096/720)e+(15625/720)g)(f'^(vi)(i) Deltax^6)+.....`

Now divide above equation by `Delta x^2` on both side

Hence we get,

`(d^2u)/dx^2 =` `(0a+0.5b+2c+(9/2)d+8e+(25/2)g)f''(i)+``(0a+(1/720)b+(64/720)c+(729/720)d+(4096/720)e+(15625/720)g)(f'^(vi)(i) Deltax^4)+.....`

Which contains `Deltax^4` i.e. 4th order and higher (above 4) order terms.

MATLAB Code for Comparison of errors:- 

clear all
close all
clc

%analytical function, f(x)= cos(x)*exp(x)
x=pi/3;
%hence, first order analytical derivative can be calculated as, f'(x)= (cos(x)*exp(x))- (sin(x)*exp(x))
%hence, second order analytical derivative can be calculated as, f''(x)= -2*(sin(x)*exp(x))
Analytical_derivative = -2*(sin(x)*exp(x));

%range for dx values
dx = linspace(pi/40,pi/400,25);

for i=1:length(dx);
  %Central differencing scheme for second order derivative at fourth order approximation,
  %formula, (a*f(x-2dx)+b*f(x-dx)+c*f(x)+d*f(x+dx)+e*f(x+2dx))/(dx)^2
  %X=[a b c d e]
  X = [-0.08333333333333337 1.333333333333334 -2.500000000000001 1.333333333333335 -0.0833333333333337]; 
  central = (X(1)*cos(x-2*dx(i))*exp(x-2*dx(i)) + X(2)*cos(x-dx(i))*exp(x-dx(i)) + X(3)*cos(x)*exp(x) + X(4)*cos(x+dx(i))*exp(x+dx(i)) + X(5)*cos(x+2*dx(i))*exp(x+2*dx(i)) )/((dx(i))^2);
  Central_differencing_error(i) = abs(Analytical_derivative - central);

  %Right skewed differencing scheme for second order derivative at fourth order approximation,
  %formula, (a*f(x)+b*f(x+dx)+c*f(x+2dx)+d*f(x+3dx)+e*f(x+4dx)+g*f(x+5dx))/(dx)^2
  %Y=[a b c d e g]
  Y = [3.749999999999991 -12.8333333333333 17.83333333333329 -12.99999999999998 5.083333333333328 -0.8333333333333329]; 
  right_skewed = (Y(1)*cos(x)*exp(x) + Y(2)*cos(x+dx(i))*exp(x+dx(i)) + Y(3)*cos(x+2*dx(i))*exp(x+2*dx(i)) + Y(4)*cos(x+3*dx(i))*exp(x+3*dx(i)) + Y(5)*cos(x+4*dx(i))*exp(x+4*dx(i)) + Y(6)*cos(x+5*dx(i))*exp(x+5*dx(i)) )/((dx(i))^2);
  Right_skewed_differencing_error(i) = abs(Analytical_derivative - right_skewed);

  %Left skewed differencing scheme for second order derivative at fourth order approximation,
  %formula, (a*f(x)+b*f(x-dx)+c*f(x-2dx)+d*f(x-3dx)+e*f(x-4dx)+g*f(x-5dx))/(dx)^2
  %Z=[a b c d e g]
  Z = [3.749999999999991 -12.8333333333333 17.83333333333329 -12.99999999999998 5.083333333333328 -0.8333333333333329]; 
  left_skewed = (Z(1)*cos(x)*exp(x) + Z(2)*cos(x-dx(i))*exp(x-dx(i)) + Z(3)*cos(x-2*dx(i))*exp(x-2*dx(i)) + Z(4)*cos(x-3*dx(i))*exp(x-3*dx(i)) + Z(5)*cos(x-4*dx(i))*exp(x-4*dx(i)) + Z(6)*cos(x-5*dx(i))*exp(x-5*dx(i)) )/((dx(i))^2);
  Left_skewed_differencing_error(i) = abs(Analytical_derivative - left_skewed);
  
endfor

figure(1);
loglog(dx,Central_differencing_error,'-d-','color','r');
hold on;
loglog(dx,Right_skewed_differencing_error,'-o-','color','b');
hold on;
loglog(dx,Left_skewed_differencing_error,'-s-','color','m');
axis([pi/450 pi/30 10^(-11) 10^(-1)]);
legend('Central differencing error','Right skewed differencing error','Left skewed differencing error');
xlabel('dx');
ylabel('Truncation error');
title('Comparison of differencing schemes');

 Given data:- Function, `f(x)=cos(x).e^x`

As per `d/dx(uxxv)= (u'xxv)+(v'xxu)`

Where `u=cos(x)` & `v=e^x`

Hence Second order derivative, `f''(x) = -2.sin(x).e^x=-4.935745353035276`

As per above derived equations, we have formulae for 

Central difference scheme, `(-0.083333*f(x-2dx) + 1.3333*f(x-dx) - 2.5*f(x) + 1.3333*f(x+dx) - 0.083333*f(x+2dx))/(dx)^2` 

Right side skewed scheme, `(3.75*f(x) - 12.833*f(x+dx) + 17.833*f(x+2dx) - 13*f(x+3dx) + 5.0833*f(x+4dx) - 0.83333*f(x+5dx))/(dx)^2`

Left side skewed scheme, `(3.75*f(x) - 12.833*f(x-dx) + 17.833*f(x-2dx) - 13*f(x-3dx) + 5.0833*f(x-4dx) - 0.83333*f(x-5dx))/(dx)^2`

So, truncation error is calculated by takin differennce between absolute derivative value and numerical approximate value.

Above plot is 'dx values' v/s 'truncation error', where plot is generated on logarithmic scale for different values of 'dx'. From above graph it is clear that the truncation error values offered by central differencing are less than right & left skewed differencing scheme i.e. Central differencing scheme gives closed solution to true solution than skewed schemes. While right side skewed error and left side skewed errors are merely equal

Hence, `"Central differencing error"  <  "Right skewed error" ~~ "Left skewed error"` 

Advantage of skewed scheme over central differencing scheme:- As skewed schemes required data from only one end i.e. either left side or right side from the considered node, so its not complicated to use skewed scheme when direction of any property is given. But if we consider central differencing scheme we have to provide data from both side of considered node, so its bit complex to provide data from both side when flow direction of property is given. Hence for such cases skewed schemes are better and more accurate than central differencing scheme.