Taylor Table method for deriving central and skewed schemes of desired accuracy

1.)Derivation of central differencing of fourth order approximation and second derivative:

Hypothesis function:

del x^2 * (d^2 f/dx^2)= a.f(i-2) + b.f(i-1) + c.f(i) + d.f(i+1) + e.f(i+2)

The \".\" above doesnt represent the elementwise operator and rather is a simple multiplication symbol used for representation in a cleaner manner.

We first go on to expand the term a.f(i-2) using taylor series and similarly all other right hand side terms subsequently.

a.f(i-2)= a[f(i) - (2del x*f\'(i)) + ((2del x)^2/2!)*f\"(i) - ((2del x)^3/3!) * f\'\'\'(i) + ((2del x)^4/4!) * f\'\'\'\'(i)]

=af(i) - 2a del x f\'(i) + 2a (del x)^2 f\'\'(i) - (8a/6)(del x)^3 f\'\'\'(i) + (16a/24) (del x)^4 f\'\'\'\'(i)

As in the above line,the terms say for eg: af(i) obviously means a*f(i) but for time saving purposes we have written it as af(i) and similarly for other terms.

bf(i-1) = b[f(i) - (del x * f\'(i)) + ((del x)^2/2!)*f\"(i) - ((del x)^3/3!)*f\'\'\'(i) + ((del x)^4/4!) * f\'\'\'\'(i)]

= bf(i) - b del x f\'(i) + b (del x^2/2) f\"(i) - b (del x^3/6) f\'\'\'(i) + b (del x^4/24) f\'\'\'\'(i)

c.f(i) = c.f(i)

d.f(i+1) = d[f(i) + (del x * f\'(i)) + (del x^2/2!)*f\'\'(i) + (del x^3 /3!)* f\'\'\'(i) + (del x^4/4!)*f\'\'\'\'(i)]

= df(i) + d delx f\'(i) + (d/2) del x^2 f\"(i) + (d/6) del x^3 f\'\'\'(i) + (d/24) del x^4 f\'\'\'\'(i)

e.f(i+2) = e[f(i) + ((2 del x)*f\'(i)) + ((2 del x)^2 / 2!)* f\'\'(i) + ((2 del x)^3/3!)*f\'\'\'(i) + ((2 del x)^4/4!)* f\'\'\'\'(i)]

= ef(i) + e 2 del x f\'(i) + e 2 del x^2 f\"(i) + (8/6) e del x ^3 f\'\'\'(i) + (16/24) e del x^4 f\'\'\'\'(i)

We then go on to form the taylor table with the numerical stencil as the 1st column and taylor series terms as the 1st row,and then go on to fill the coefficients for the respective terms of taylor series of each expansion done. Since we have five unknowns,we consider five equations only to solve for a , b , c , d , e. The equations formed are :

a+b+c+d+e = 0

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

2a+b/2 + d/2 + 2e=1

-4/3a -b/6 + d/6 + 4/3e = 0

2/3a + b/24 + d/24 + 2/3e=0

Since, we are computing the second order serivative hence the sum of the coefficients of the column del x^2 f\"(i) = 1 and all the rest=0.

The matrix form is obtained from the above five equations and it would look like:

[A][X]= [B]

where [A]= [1 1 1 1 1;-2 -1 0 1 2; 2 1/2 0 1/2 2; -4/3 -1/6 0 1/6 4/3; 2/3 1/24 0 1/24 2/3] where \";\" separates each row from one another , [X]= [a;b;c;d;e] and [B]= [0;0;1;0;0]. To find the solution [X],we do [A]^-1 * [B] where [A]^-1 can be called as A inverse and is computed in matlab using an in built function \"inv\" and is used as inv(A).

After computation of solution:

a=-0.0833,b=1.3333,c=-2.5000,d=1.3333,e=-0.0833

Substituing the above value of coefficients in main hypothesis function we get:

(d^2 f/dx^2) = (-0.0833f(i-2) + 1.3333 f(i-1) - 2.5000f(i) + 1.3333f(i+1) - 0.0833f(i+2))/ (del x ^ 2) .

Coding :

We have created a function named as \"central_diff\"that can take in x,dx values as inputs and output an error term which is assigned to a variable \"op_challenge\". The analytical function which we have considered is: (exp (x) * cos(x)). The second derivative of the analytical function calculated manually is given by : - 2 *sin(x) *exp(x) is stored in\"analytical_derivative_second\". The matrix form is then entered represented by [A],[B] and [X] where [X] ouputs the solution for unknown coefficients a,b,c,d,e and consists of 5 rows and 1column. The formula that has been computed using taylor series is then keyed in into the function using the variable central_scheme and the absolute difference between the central_scheme and analytical_derivative_second is stored in the output of the function which is\"op_challenge\". The \"function\" keyword begins the function and the end keyword ends it.

function [op_challenge] = central_diff(x,dx)
%Analytical function
%f(x)= (exp(x)*cos(x))
%Computing second order derivative of analytical function
analytical_derivative_second = -2 * sin(x)*exp(x)
A= [1 1 1 1 1 ; -2 -1 0 1 2; 2 1/2 0 1/2 2; -4/3 -1/6 0 1/6 4/3;  2/3 1/24 0 1/24 2/3]
B=[0;0;1;0;0]
X= inv(A)*B;

%Computation of numerical derivative
%central_scheme_fourth_order = (-0.0833*f(x-(2*dx)) + 1.3333*f(x-dx)
%-2.5000*f(x) + 1.3333 *f(x+dx) - 0.0833*f(x+(2*dx)))/(dx^2);

central_scheme = ((X(1) * exp(x-(2*dx)) *cos(x-(2*dx))) + (X(2)*exp(x-dx)*cos(x-dx)) + (X(3) * exp(x)* cos(x)) + (X(4) * exp(x+dx)*cos(x+dx)) + (X(5)*exp(x+(2*dx))*cos(x+(2*dx))))/(dx^2)

op_challenge=abs(central_scheme - analytical_derivative_second)
end

2.)Derivation of skewed right side differencing of fourth order approximation and second derivative:

Hypothesis function:

del x^2 * (d^2 f/d ^2)= a.f(i) + b.f(i+1)+ c.f(i+2) +d.f(i+3) + e.f(i+4) + g.f(i+5)

The \".\" above doesnt represent the elementwise operator and rather is a simple multiplication symbol used for representation in a cleaner manner.

We first go on to expand the term b.f(i+1) and all other right hand side terms subsequently.

b.f(i+1)= b[f(i) + (del x* f\'(i)) + (del x^2/2!)*f\"(i) + (del x^3/3!)* f\'\'\'(i) + (del x^4/4!)*f\'\'\'\'(i) + (del x^5/5!)*f\'\'\'\'\'(i)]

= bf(i) + b del x f\'(i) + b/2 del x^2 f\"(i) + b/6 del x^3 f\'\'\'(i) + b/24 del x^4 f\'\'\'\'(i) + b/120 del x^5 f\'\'\'\'\'(i)

As in the above line,the terms say for eg: b del x f\'(i) obviously means b*del x* f\'(i) but for time saving purposes we have written it as b del x f\'(i) and similarly for other terms.

c.f(i+2) = c[f(i) + (2 del x)* f\'(i) + ((2 del x)^2/2!) * f\'\'(i) + ((2 del x)^3/3!)*f\'\'\'(i) +((2 del x)^4/4!)* f\'\'\'\'(i) + ((2 del x)^5 /5!)* f\'\'\'\'\'(i)]

= cf(i)+c (2 del x f\'(i)) + c(2 del x^2 f\'(i)) + c(8/6 del x^3 f\'\'\'(i)) + c(16/24 del x^4 f\'\'\'\'(i)) + c(32/120 del x^5 f\'\'\'\'\'(i))

d.f(i+3) = d[f(i) + (3 del x * f\'(i)) + ((3 del x)^2/2!)*f\"(i) + ((3 del x)^3/3!)*f\'\'\'(i) + ((3 del x)^4/4!)*f\'\'\'\'(i) + ((3 del x)^5/5!)*f\'\'\'\'\'(i)]

= df(i) + d (3 del x f\'(i)) + d(9/2 del x^2 f\"(i)) + d(27/6 del x^3 f\'\'\'(i)) + d(81/24 del x^4 f\'\'\'\'(i)) + d( 243/120 del x^5 f\'\'\'\'\'(i))

e.f(i+4) = e[f(i) + (4 del x)*f\'(i) + ((4 del x)^2/2!)*f\"(i) + ((4 del x)^3/3!) *f\'\'\'(i) + ((4 del x)^4/4!)*f\'\'\'\'(i) + ((4 del x)^5/5!)*f\'\'\'\'\'(i)]

= ef(i) + e(4 del x f\'(i)) + e(8 del x^2 f\'\'(i)) + e(32/3 del x^3 f\'\'\'(i)) + e(32/3 del x^4 f\'\'\'\'(i)) + e(128/15 del x^5 f\'\'\'\'\'(i))

g.f(i+5) = g[f(i) + (5 del x)*f\'(i) + ((5 del x)^2/2!)*f\"(i) + ((5 del x)^3/3!)*f\'\'\'(i) + ((5 del x)^4/4!)* f\'\'\'\'(i) + ((5 del x)^5/5!)*f\'\'\'\'\'(i)]

=gf(i) + g(5 del x f\'(i)) + g(25/2 del x^2 f\"(i)) + g(125/6 del x^3 f\'\'\'(i)) + g(625/24 del x^4 f\'\'\'\'(i)) + g(3125/120 del x^5 f\'\'\'\'\'(i))

a.f(i) = a.f(i)

Also,note that the term for eg: g(125/6 del x^3 f\'\'\'(i)) means g*(125/6)*(del x^4)*f\'\'\'(i) and similarly for other terms .

We then go on to form the taylor table with the numerical stencil as the 1st column and taylor series terms as the 1st row,and then go on to fill the coefficients for the respective terms of taylor series of each expansion done. Here,we basically have six unknowns and hence six equations in total which are given below:

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

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

b/2 + 2c + (9/2)d + 8e (25/2)g=1

b/6 + (4/3)c + (9/2)d + (32/3)e + (125/6)g = 0

b/24 +(2/3)c + (27/8)d + (32/3)e + (625/24)g = 0

b/120 + (4/15)c + (243/120)d + (128/15)e + (625/24)g=0

Since, we are computing the second order serivative hence the sum of the coefficients of the column del x^2 f\"(i) = 1 and all the rest=0.

The matrix form is obtained from the above six equations and it would look like:

[A][X]= [B]

[A]= [1 1 1 1 1 1;0 1 2 3 4 5;0 1/2 2 9/2 8 25/2;0 1/6 4/3 9/2 32/3 125/6; 0 1/24 2/3 27/8 32/3 625/24; 0 1/120 4/15 243/120 128/15 625/24] where \";\" separates each row from one another , [X]= [a;b;c;d;e;g] and [B]= [0;0;1;0;0;0]. To find the solution [X],we do [A]^-1 * [B] where [A]^-1 can be called as A inverse and is computed in matlab using an in built function \"inv\" and is used as inv(A).

After computation of solution:

a=-3.7500,b=-12.8333,c=17.8333,d=-13.0000,e=5.0833,g= -0.8333

Substituing the above value of coefficients in main hypothesis function we get:

(d^2 f/dx^2) = (-3.7500f(i) - 12.8333 f(i+1) + 17.8333f(i+2) - 13.0000f(i+3) + 5.0833f(i+4) - 0.8333f(i+5))/ (del x ^ 2) .

Coding :

We have created a function named as \"skewed_rightside\"that can take in x,dx values as inputs and output an error term which is assigned to a variable \"op_challenge\". The analytical function which we have considered is: (exp (x) * cos(x)). The second derivative of the analytical function calculated manually is given by : - 2 *sin(x) *exp(x) is stored in\"analytical_derivative_second\". The matrix form is then entered represented by [A1],[B1] and [X1] where [X1] ouputs the solution for unknown coefficients a,b,c,d,e,g and consists of 6 rows and 1column. The formula that has been computed using taylor series is then keyed in into the function using the variable skewedright_scheme and the absolute difference between the skewedright_scheme and analytical_derivative_second is stored in the output of the function which is\"op_challenge\".

function [op_challenge] = skewed_rightside(x,dx)
%Analytical function
%f(x)= (exp(x)*cos(x))
%Computing second order derivative of analytical function
analytical_derivative_second = -2 * sin(x)*exp(x)
A1= [1 1 1 1 1 1; 0 1 2 3 4 5; 0 1/2 2 9/2 8 25/2; 0 1/6 4/3 9/2 32/3 125/6; 0 1/24 2/3 27/8 32/3 625/24; 0 1/120 4/15 243/120 128/15 625/24]
B1=[0;0;1;0;0;0]
X1= inv(A1)*B1;

%Computation of numerical derivative
%skewed_right_sided_scheme = (-3.7500*f(x)) - (12.8333*f(x+dx)) +
%((17.8333*f(x+(2*dx))) - (13.0000*f(x+(3*dx))) + (5.0833*f(x+(4*dx))) -
%(0.8333*f(x+(5*dx)))))/(dx^2)

skewedright_scheme = ((X1(1)*exp(x)*cos(x)) +(X1(2)*exp(x+dx)*cos(x+dx)) + (X1(3)*exp(x+(2*dx))*cos(x+(2*dx))) +(X1(4) *exp(x+(3*dx))*cos(x+(3*dx))) + (X1(5)*exp(x+(4*dx))*cos(x+(4*dx)))+(X1(6)*exp(x+(5*dx))*cos(x+(5*dx))))/(dx^2)
op_challenge=abs(skewedright_scheme - analytical_derivative_second)
end

3.)Derivation of skewed left side differencing of fourth order approximation and second derivative:

Hypothesis function:

del x^2 * (d^2 f/d ^2)= a.f(i-5) + b.f(i-4)+ c.f(i-3) +d.f(i-2) + e.f(i-1) + g.f(i)

The \".\" above doesnt represent the elementwise operator and rather is a simple multiplication symbol used for representation in a cleaner manner.

We first go on to expand the term a.f(i-5) and all other right hand terms subsequently.

a.f(i-5) = a[f9i) - (5 del x)*f\'(i) + ((5 del x)^2/2!)*f\"(i) - ((5 del x)^3/3!)*f\'\'\'(i) + ((5 del x)^4/4!) * f\'\'\'\'(i) - ((5 del x)^5/5!)*f\'\'\'\'\'(i)]

= af(i) - a(5 del x f\'(i)) + a(25/2 del x^2 f\"(i)) - a(125/6 del x^3 f\'\'\'(i)) + a(625/24 del x^4 f\'\'\'\'(i)) - a(3125/120 del x^5 f\'\'\'\'\'(i))

As in the above line,the terms say for eg: -a(25/2 del x ^2 f\'\'(i)) obviously means - a*del x^2* f\'\'(i)* (25/2) but for time saving purposes we have written it as - a(25/2 del x^2 f\'\'(i)) and similarly for other terms.

b.f(i-4) = b[f(i) - (4 del x)*f\'(i) + ((4 del x)^2/2!)*f\'\'(i) - ((4 del x)^3/3!) * f\'\'\'(i) + ((4 del x)^4/4!)*f\'\'\'\'(i) - ((4 del x)^5/5!)*f\'\'\'\'\'(i)]

= bf(i) - b(4 del x f(i)) + b(8 del x^2 f\"(i)) - b(32/3 del x^2 f\'\'\'(i)) + b(32/3 del x^4 f\'\'\'\'(i)) - b(128/15 del x^5 f\'\'\'\'\'(i))

c.f(i-3) = c[f(i) - (3 del x)*f\'(i) + ((3 del x)^2/2!)* f\"(i) - ((3 del x)^3/3!) * f\'\'\'(i) + ((3 del x)^4/4!)*f\'\'\'\'(i) - ((3 del x)^5/5!)*f\'\'\'\'\'(i)]

= cf(i) - c(3 del x f\'(i)) + c(9/2 del x^2 f\"(i)) - c(9/2 del x^3 f\'\'\'(i)) + c(27/8 del x^4 f\'\'\'\'(i)) - c(243/120 del x^5 f\'\'\'\'\'(i))

d.f(i-2) = d[f(i) - (2 del x)* f\'(i) + ((2 del x)^2/2!)* f\'\'(i) - ((2 del x)^3/3!)*f\'\'\'(i) + ((2 del x)^4/4!)* f\'\'\'\'(i) - ((2 del x)^5/5!)* f\'\'\'\'\'(i)]

= df(i) - d(2 del x f\'(i)) + d(2 del x^2 f\"(i)) - d(4/3 del x^3 f\'\'\'(i)) + d(2/3 del x^4 f\'\'\'\'(i)) - d(32/120 del x^5 f\'\'\'\'\'(i))

e.f(i-1) = e[f(i) - (del x * f\'(i)) + ((del x)^2/2!)*f\"(i) - ((del x)^3/3!)* f\'\'\'(i) + ((del x)^4/4!)*f\'\'\'\'(i) - ((del x)^5/5!)*f\'\'\'\'\'(i)]

= ef(i) - e(del x f\'(i)) + e(1/2 del x^2 f\'\'(i)) - e(1/6 del x^3 f\'\'\'(i)) + e(1/24 del x^4 f\'\'\'\'(i)) - e(1/120 del x^5 f\'\'\'\'\'(i))

g.f(i) = g.f(i)

We then go on to form the taylor table with the numerical stencil as the 1st column and taylor series terms as the 1st row,and then go on to fill the coefficients for the respective terms of taylor series of each expansion done. Here,we basically have six unknowns and hence six equations in total which are given below:

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

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

(25/2)a + 8b +(9/2)c + 2d + (1/2)e=1

(-125/6) a - (32/3)b - (9/2)c -(4/3)d - (1/6)e= 0

(625/24)a + (32/3)b + (27/8)c + (2/3)d + (1/24)e = 0

(-625/24)a -(128/15)b -(243/120)c - (4/15)d - (1/120)e = 0

Since, we are computing the second order serivative hence the sum of the coefficients of the column del x^2 f\"(i) = 1 and all the rest=0.

The matrix form is obtained from the above six equations and it would look like:

[A][X]= [B]

[A]= [1 1 1 1 1 1;-5 -4 -3 -2 -1 0; 25/2 8 9/2 2 1/2 0;-125/6 -32/3 -9/2 -4/3 -1/6 0;625/24 32/3 27/8 2/3 1/24 0; -625/24 -128/15 -243/120 -4/15 1/120 0] where \";\" separates each row from one another , [X]= [a;b;c;d;e;g] and [B]= [0;0;1;0;0;0]. To find the solution [X],we do [A]^-1 * [B] where [A]^-1 can be called as A inverse and is computed in matlab using an in built function \"inv\" and is used as inv(A).

After computation of solution:

a=-1.0667,b=6.2500,c=-15.3333,d=20.1667,e=-14.0000,g= 3.9833

Substituing the above value of coefficients in main hypothesis function we get:

(d^2 f/dx^2) = (-1.0667f(i-5) + 6.2500 f(i-4) -15.3333f(i-3) + 20.1667f(i-2) - 14.0000f(i-1) + 3.9833f(i))/ (del x ^ 2)

Coding :

We have created a function named as \"skewed_leftside\"that can take in x,dx values as inputs and output an error term which is assigned to a variable \"op_challenge\". The analytical function which we have considered is: (exp (x) * cos(x)). The second derivative of the analytical function calculated manually is given by : - 2 *sin(x) *exp(x) is stored in\"analytical_derivative_second\". The matrix form is then entered represented by [A2],[B2] and [X2] where [X2] ouputs the solution for unknown coefficients a,b,c,d,e,g and consists of 6 rows and 1column. The formula that has been computed using taylor series is then keyed in into the function using the variable skewedleft_scheme and the absolute difference between the skewedleft_scheme and analytical_derivative_second is stored in the output of the function which is\"op_challenge\".

function [op_challenge] = skewed_leftside(x,dx)
%Analytical function
%f(x)= (exp(x)*cos(x))
%Computing second order derivative of analytical function
analytical_derivative_second = -2 * sin(x)*exp(x)
A2=[1 1 1 1 1 1; -5 -4 -3 -2 -1 0; 25/2 8 9/2 2 1/2 0; -125/6 -32/3 -9/2 -4/3 -1/6 0; 625/24 32/3 27/8 2/3 1/24 0; -625/24 -128/15 -243/120 -4/15 1/120 0]
B2=[0;0;1;0;0;0]
X2= inv(A2)* B2

%Computation of numerical derivative
%skewed_left_sided_scheme = ((-1.0667*f(x-(5*dx)))+ (6.2500*f(x-(4*dx))) -
%(15.3333*f(x-(3*dx))) + (20.1667*f(x-(2*dx))) - (14.0000*f(x-dx)) +
%(3.9833*f(x)))/(dx^2)

skewedleft_scheme = ((X2(1)*exp(x-(5*dx))*cos(x-(5*dx))) + (X2(2)*exp(x-(4*dx))*cos(x-(4*dx))) + (X2(3)*exp(x-(3*dx))*cos(x-(3*dx))) + (X2(4)*exp(x-(2*dx))*cos(x-(2*dx))) +(X2(5)*exp(x-dx)*cos(x-dx)) + (X2(6)*exp(x)*cos(x)))/(dx^2)
op_challenge=abs(skewedleft_scheme - analytical_derivative_second)
end

3 (a)Plot description:

As we can see from the above plot that the error is the least in the central differencing scheme as compared to the skewed right side and skewed left side scheme ,for the same order of accuracy i.e fourth order of the second order derivative,we can say that the plot is correct. The \'blue\" line represents the absolute error value between second derivative computed analytically and central differencing numerical scheme and how it changes for a range of dx. The \'red\' and \'green\' line shows the absolute error values computed through skewed right side and skewed left side schemes respectively.

We have also noticed that for a range of dx value upto say 10^-2,the behaviour of the various order approximations turn out to be of the trend :

central_error<skewedr_error~skewedl_error which is correct. But,as we make dx even smaller, the error starts to increase for all the schemes due to the increase in complication in calculation. Hence, the right grid spacing \'dx\'and range of dx would be very essential while deriving accurate schemes with utmost efficiency.

Coding last part:

Here, we have considered x=pi/3 and 100 different values of dx between pi/4 to pi/4000 using dx= linspace(pi/4,pi/4000,100). We have then run a \'for loop\' upto the 100th value of dx which is done using length(dx). i represents the index or position of each value . For eg: central_error(i) in the first loop will consider the first value of dx and compute central_error,and similarly for skewedr_error and skewedl_error.

The plotting has been done using \'loglog\' command due to the disparity in magnitudes and gives a much better representation as compared to \'plot\' command.

The legends have been provided using legend command and xlabel,y label as well.

clear all
close all
clc

x= pi/3
dx= linspace(pi/4,pi/4000,100)

for i=1:length(dx)
    
    central_error(i)=central_diff(x,dx(i))
    skewedr_error(i)= skewed_rightside(x,dx(i))
    skewedl_error(i) = skewed_leftside(x,dx(i))
end

figure(1)
loglog(dx, central_error)
hold on

loglog(dx, skewedr_error,\'r\')
loglog(dx, skewedl_error,\'g\')
xlabel(\'Range of dx\')
ylabel(\'Absolute Errors\')
legend(\'Central differencing of fourth prder approx\',\'Right sided skewed fourth order approx\',\'Left sided Skewed fourth order approx\')

Note:

We exactly knew how many terms to consider in the hypothesis function due to the help from challenge link where it is said : No of nodes = p+q where p is the order of approximation which in this case is 4,and and q= order of derivative which is 2 . This holds true for skewed schemes only . Hence,we consider 6 terms exactly. For Central differencing schemes, the formula is p+q - 1 with the same convention. Hence,the no of terms we consider is 5.

(3b) A skewed scheme can generate the required order of accuracy for a non- symmetric stencil and includes 1 more point as compared to central scheme whreas a central scheme is limited to symmetric type of stencils only.