Week 3 - Taylor table method and Matlab code

Derivation for 4th order approximations of the second-derivative:

1. Central difference:

For the central differencing, there will be a stencil group consisting of 5 nodal points of the ith node in order to derive the 4th order approximations of the 2nd derivative as follows

so the R.H.s of the equation of the second derivative will be `(d^2u) / (dx^2)`

for the central difference scheme,

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

where a,b,c,d,e are constants and are to be calculated in this derivation.

A Taylor table is being used for the calculation as below

	f(i)	dxf'(i)	`dx^2f''(x)`	`dx^3f'''(x)`	`dx^4f^(iv)(x)`	`dx^5f^(v)(x)`	`dx^6f^(vi)(x)`
af(i-2)	a	-2a	4a/2	-8a/6	16a/24	-32a/120	64a/740
bf(i-1)	b	-b	b/2	-b/6	b/24	-b/120	b/740
cf(i)	c	0	0	0	0	0	0
df(i+1)	d	d	d/2	d/6	d/24	d/120	d/740
ef(i+2)	e	2e	4e/2	8e/6	16e/24	32e/120	64e/740
`sum`	0	0	1	0	0	?	?

The table has been arranged according to the Taylor series expansion and the sum which will be equal to R.H.S i.e., `(d^2u) / (dx^2)`

As the sum has to be equivalent to R.H.S, all the second derivative terms i.e., f''(x) will be equivalent to one and the ones that don't should be zero.

Considering the 5th and 6th order derivative terms the value is unknown as we have the equation for 5 terms. specifying the 6th and 7th as zero would simply be an over-simplified system.

summing the terms, we get the following equations

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

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

`4a/2 +b/2 + 0 + d/2 + 4e/2 = 1`

`-8a/6 -b/6 + 0 +d/6+ 8e/6=0`

`16a/24+b/24+0+d/24+16e/24=0`

now putting in matrix form,

`A = [[1,1,1,1,1],[-2,-1,0,1,2],[2,1/2,0,1/2,2],[-8/6,-1/6,0,1/6,8/6],[16/24,1/24,0,1/24,16/24]]`

`X= [[a], [b], [c], [d],[e]]`

`B= [[0],[ 0],[ 1],[ 0],[ 0]]`

where A. X = B;

using inverse square law,

X= inv(A) * B

`[[a], [b], [c], [d],[e]] = [[-0.0833],[1.3333], [-2.5000],[1.3333],[-0.0833]]`

by comparing the values, we could say that, [a=e] and [b=d]

now simplifying the terms, the equation we get

af(i-2) + bf(i-1)+cf(i)+df(i+1)+ef(i+2)` = `(2a+b/2+0+d/2+2e)f''(x) dx^2 + (64a/240 + b/240+0+d/240+64e/240)f^(1v) dx^6`

`af(i-2) + bf(i-1)+cf(i)+df(i+1)+ef(i+2)` = `(2a+b/2+0+d/2+2e)f'' dx^2 + (0) f^(v1) dx^6`

dividing by `dx^2`, we get

`(d^2u)/(dx^2)``~=` `(af(i-2) + bf(i-1)+cf(i)+df(i+1)+ef(i+2))/dx^2` = `(2a+b/2+0+d/2+2e)f'' + (0) f^(v1) dx^4`

Therefore, This scheme is central difference for 4th order approximation for second order derivative.

2. Skewed right-sided difference:

For the Skewed right-sided difference, there will be a stencil group consisting of 5 nodal points of the ith node at the start point to the end from the right side, the following sequence of the ith order will be as shown below:

so the R.H.s of the equation of the second derivative will be equivalent to `(d^2u) / (dx^2)`

for Skewed right-sided difference,

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

where a,b,c,d,e are constants and are to be calculated in this derivation.

The following Taylor expansion is done using Taylor table as below

	f(i)	dxf'(i)	`dx^2f''(x)`	`dx^3f'''(x)`	`dx^4f^(iv)(x)`	`dx^5f^(v)(x)`	`dx^6f^(vi)(x)`
`af(i)`	a	0	0	0	0	0	0
`bf(i+1)`	b	b	b/2	b/6	b/24	b/120	b/740
`cf(i+2)`	c	2c	4c/2	8c/6	16c/24	32c/120	64c/740
`df(i+3)`	d	3d	9d/2	27d/6	81d/24	243d/120	729d/740
`ef(i+4)`	e	4e	16e/2	64e/6	256e/24	1024e/120	4096e/740
gf(i+5)	g	5g	25g/2	125g/6	625g/24	3125g/120	15625g/740
`sum`	0	0	1	0	0	0	?

As the 7th terms will be unknown and the summing the terms, we get the following equations for constants;

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

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

`0+b/2+4c/2+9d/2+16e/2+25g/2=1`

`0+b/6+8c/6+27d/6+64e/6+125g/6=0`

`0+b/24+16c/24+81d/24+256e/24+625g/24=0`

`0+b/120+32c/120+243d/120+1024e/120+3125g/120 = 0`

Arranging these in matrix form, we get

`A= [[1,1,1,1,1,1],[0,1,2,3,4,5],[0,1/2,4/2,9/2,16/2,25/2],[0,1/6,8/6,27/6,64/6,125/6],[0,1/24,16/24,81/24,256/24,625/24],[0,1/120,32/120,243/120,1024/120,3125/120]]`

`X= [[a], [b], [c], [d],[e],[g]]`

`B= [[0],[ 0],[ 1],[ 0],[ 0],[0]]`

where A. X = B;

using inverse square law in matlab, we calculate the X

X = inv(A)* B

`[[a], [b], [c], [d],[e],[g]] = [[3.7500],[-12.8333], [17.3333],[-13.0000],[5.0833],[-0.8333]]`

Now, simplifying the terms with constants, we get

`af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(1+4)+gf(i+5)` = `(0+b/2+4c/2+9d/2+16e/2+25g/2)dx^2 f''(x)+ (0+b/740+64c/740+729d/740+4096e/740+15625g/740)dx^6 f^(vi)(x)`

dividing with dx^2, we get

`(d^2u) / (dx^2)``~=``(af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(1+4)+gf(i+5))/dx^2`=`(0+b/2+4c/2+9d/2+16e/2+25g/2) f''(x)+ (0)dx^4 f^(vi)(x)`

which is., this skewed right sided scheme is 4th order approximation for second order derivative.

3. Skewed left-sided difference:

For the Skewed Left-sided difference, there will be a stencil group consisting of 5 nodal points of the ith node at the start point to the end from the left side, the following sequence of the ith order will be as shown below:

so the R.H.s of the equation of the second derivative will be equivalent to `(d^2u) / (dx^2)`

for Skewed left-sided difference,

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

where a,b,c,d,e are constants and are to be calculated in this derivation.

The following Taylor expansion is done using Taylor table as below

	f(i)	`dx^2f'(x)`	`dx^2f''(x)`	`dx^3f'''(x)`	`dx^4f^(iv)(x)`	`dx^5f^(v)(x)`	`dx^6f^(vi)(x)`
af(i-5)	a	-5a	25a/2	-125a/6	652a/24	-3125a/120	15625a/740
bf(i-4)	b	-4b	16b/2	-64b/6	256b/24	-1024b/120	4096b/740
cf(i-3)	c	-3c	9c/2	-27c/6	81c/24	-243c/120	729c/740
df(i-2)	d	-2d	4d/2	-8d/6	16d/24	-32d/120	64d/740
ef(1-1)	e	-e	e/2	-e/6	e/24	-e/120	e/740
gf(i)	g	0	0	0	0	0	0
`sum`	0	0	1	0	0	0	?

As the 7th terms will be unknown and specifying it will lead to over-specification of the system

Summing the terms, we get the following equations for constants;

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

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

`25a/2+16b/2+9c/2+4d/2+e/2+0=1`

`-125a/6-64b/6-27c/6-8d/6-e/6-0=0`

`652a/24+256b/24+81c/24+16d/24+e/24+0=0`

`-3125a/120-1024b/120-243c/120-32d/120-e/120+0 = 0`

Arranging these equations of constant in matrix form, we get

`A = [ [1,1,1,1,1,1], [-5,-4,-3,-2,-1,0], [25/2,16/2,9/2,4/2,1/2,0], [-125/6,-64/6,-27/6,-8/6,-1/6,0],[652/24,256/24,81/24,16/24,1/24,0],[-3125/120,-1024/120,-243/120,-32/120,-1/120,0]]`

`X= [[a], [b], [c], [d],[e],[g]]`

`B= [[0],[ 0],[ 1],[ 0],[ 0],[0]]`

using inverse square law in matlab, we calculate the X

X = inv(A)* B

Therefore,`[[a], [b], [c], [d],[e],[g]] = [[-0.8333],[5.0833],[-13.000],[17.8333],[-12.8333],[3.7500]]`

Now, simplifying the terms with constants, we get

`af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(1-1)+gf(i)`= `(25a/2+16b/2+9c/2+4d/2+e/2+0)dx^2f''(x) +(15625a/740+4096b/740+729c/740+64d/740+e/740+0)dx^6f^(vi)(x)`

dividing with dx^2,we get

`(d^2u) / (dx^2) ~= (af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(1-1)+gf(i))/dx^2 = (25a/2+16b/2+9c/2+4d/2+e/2+0)f''(x) +(0)dx^4f^(vi)(x)`

Therefore, This skewed left sided scheme is 4th order approximation for second order derivative.

since, the constants of both Skewed right-sided difference and Skewed left-sided difference are same but in reverse order, we could say these schemes are fourth-order accurate.

Coding a plot that compares the absolute error and the different derived schemes:

In this code, there is a plot for error and different schemes which are central differencing, skewed right-sided scheme, skewed left-sided scheme, forward and backward differencing scheme for second-order derivative as below:

clear all;
close all;
clc

x= 2*pi;
dx= linspace(pi/4,pi/400,20);

for i = 1:length(dx)
    Central_differencing_error(i)  = central_differencing(x,dx(i));
    skewed_right_sided_error(i) = skewed_right_sided(x,dx(i));
   skewed_left_sided_error(i) = skewed_left_sided(x,dx(i));
   forward_differencing_error(i)= forward_differencing(x,dx(i));
   backward_differencing_error(i)= backward_differencing(x,dx(i));
end

figure(1)
loglog(dx,Central_differencing_error)
hold on 
loglog(dx,skewed_right_sided_error,'r')
hold on
loglog(dx,skewed_left_sided_error,'y')
hold on
loglog(dx,forward_differencing_error,'g')
hold on
loglog(dx,backward_differencing_error,'m')
legend('central difference error','skewed right sided error','skewed left sided error','forward_differencing_error','backward_differencing_error')
xlabel('dx')
ylabel('absolute Error')
title('error vs different schemes')

code for function of central difference scheme:

function out = central_differencing(x,dx)

%f''(x)=-2*sin(x)*exp(x)
analytical_derivative = -2*sin(x)*exp(x);
a=-0.0833;
b= 1.3333;
c=-2.5000;
d= 1.3333;
e=-0.0833;

central_differencing = (a*fun_x(x-2*dx) + b*fun_x(x-dx)+c*fun_x(x)+d*fun_x(x+dx)+e*fun_x(x+2*dx))/dx^2;



out = abs(central_differencing -  analytical_derivative);


end

code for the function of skewed right sided scheme:

function out = skewed_right_sided(x,dx)

%f''(x)=-2*sin(x)*exp(x)
analytical_derivative = -2*sin(x)*exp(x);
a=  3.7500;
b=  -12.8333;
c=  17.3333;
d=  -13.0000;
e=  5.0833;
g=  -0.8333;

skewed_right_sided = (a*fun_x(x)+b*fun_x(x+dx)+c*fun_x(x+2*dx)+d*fun_x(x+3*dx)+e*fun_x(x+4*dx)+g*fun_x(x+5*dx))/dx^2;



out = abs(skewed_right_sided -  analytical_derivative);


end

code for the function of the skewed left-sided scheme:

function out = skewed_left_sided(x,dx)

%f''(x)=-2*sin(x)*exp(x)
analytical_derivative = -2*sin(x)*exp(x);
a=  -0.8333;
b=  5.0833;
c=  -13.0000;
d=  17.3333;
e=  -12.8333;
g=  3.7500;

skewed_left_sided = (a*fun_x(x-5*dx)+b*fun_x(x-4*dx)+c*fun_x(x-3*dx)+d*fun_x(x-2*dx)+e*fun_x(x-dx)+g*fun_x(x))/dx^2;



out = abs(skewed_left_sided -  analytical_derivative);


end

code for the function of forward differencing:

function out = forward_differencing(x,dx)

%f''(x)=-2*sin(x)*exp(x)
analytical_derivative = -2*sin(x)*exp(x);

forward_differencing = (fun_x(x+2*dx)- 2*fun_x(x+dx)+ fun_x(x))/dx^2;



out = abs(forward_differencing -  analytical_derivative);


end

code for the function of backward differencing:

function out = backward_differencing(x,dx)

%f''(x)=-2*sin(x)*exp(x)
analytical_derivative = -2*sin(x)*exp(x);

backward_differencing = (fun_x(x-2*dx)- 2*fun_x(x-dx)+ fun_x(x))/dx^2;



out = abs(backward_differencing -  analytical_derivative);


end

the function of the derivative

function out =fun_x(x)
out=exp(x)*cos(x);
end

running the following code in Matlab, the following output for the plot is:

from this plot, we see that the error for the central differencing scheme of the second-order derivative is lower compare to other schemes and the using this scheme for solving PDE's will be highly accurate.

we also see that the forward and backward differencing of the 1st order derivative have the low error compare to right and left-sided skewed schemes due to a possibility of the round off error applied in the code and also it is directly proportional to the increase in the order of error. therefore, using high order schemes always is not preferred.

The importance of the skewed schemes is that when a function needs the point of information are not available, skewed schemes comes into place to take the points of information to the farther points as normal differencing schemes takes the points of information only from the neighbours. so skewed schemes should be used and important in some cases.