Week 3 - Taylor table method and Matlab code

1. Problem definition - In this particular problem we have to first derive the 4th order approximation of second derivative using central differencing,skewed left and right sided differencing and then compute it in octave for a function exp(x)*cos(x) and analyse the error results using the plot for the above three approximations.

2. Procedure - For deriving the 4th order approximation we will make use of taylors series expansion and then we will make a taylors table for the above mentioned scheme and then we will find the coefficients of the function at various node so that we get the final approximation formula and then we will write the code in octave to compute the error value for different schemes.For writing the code we will use function calling command in octave

3. Derivation -

• CENTRAL DIFFERENCING SCHEME

`frac{del^2u}{delx^2}=frac{af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)}{Deltax^2}`

drawing taylors table-

	`f(i)`	`Deltaxf'(i)`	`Deltax^2f''(i)`	`Deltax^3f'''(i)`	`Deltax^4f''''(i)`
`af(i-2)`	a	-2a	4a/2	-8a/6	16a/24
`bf(i-1)`	b	-b	b/2	-b/6	b/24
`cf(i)`	c	0	0	0	0
`df(i+1)`	d	d	d/2	d/6	d/24
`ef(i+2)`	e	2e	2e	8e/6	16e/24
sum of all the coefficients column vise `frac{del^2u}{delx^2}`	0	0	1	0	0

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

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

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

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

`frac{16a}{24}+b/24+0+d/24+frac{16e}{24}=0`

`[[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]]*[[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`

AX = B ........(6)

where A = coefficient matrix

X = unknown matrix

B = constant matrix

premultiplying by `A^-1`

`A^-1 AX = A^-1 B`

`IX=A^-1B`

`X=A^-1B`

multiplying the matrix B with inverse of matrix A

`[[a],[b],[c],[d],[e]]=[[0,0.08333,-0.08333,-0.5,1],[0,-0.66667,1.33333,1,-4],[1,0,-2.5,0,6],[0,0.66667,1.33333,-1,-4],[0,-0.08333,-0.08333,0.5,1]]*[[0],[0],[1],[0],[0]]`

on solving the above matrix multiplication we get coefficients values as

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

therefore our final formula for second derivative 4th order approximation is -

`frac{del^2u}{delx^2}=frac{-0.0833f(i-2)+1.3333f(i-1)-2.5f(i)+1.3333f(i+1)-0.0833f(i+2)}{Deltax^2}`

The approximation is called 4th order because the order of the truncation error is 4

• SKEWED RIGHT SIDE 4th ORDER APPROXIMATION

`frac{del^2u}{delx^2}=frac{af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)}{Deltax^2}`

drawing taylors table-

	`f(i)`	`Deltaxf'(i)`	`Deltax^2f''(i)`	`Deltax^3f'''(i)`	`Deltax^4f''''(i)`	`Deltax^5f''''(i)`	`Deltax^6f'''''(i)`
`af(i)`	a	0	0	0	0	0	0
`bf(i+1)`	b	b	b/2	b/6	b/24	b/120	b/720
`cf(i+2)`	c	2c	2c	8c/6	16c/24	32c/120	64c/720
`df(i+3)`	d	3d	9d/2	27d/6	81d/24	243d/120	729d/720
`ef(i+4)`	e	4e	8e	64e/6	256e/24	1024e/120	4096e/720
`gf(i+5)`	g	5g	25g/2	125g/6	625g/24	3125g/120	15625g/720
sum of all the coefficients column vise `frac{del^2u}{delx^2}`	0	0	1	0	0	0	?

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

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

`0+b/2+2c+frac{9d}{2}+8e+frac{25g}{2}=1`

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

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

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

`0+b/720+frac{64c}{720}+frac{729d}{720}+frac{4096e}{720}+frac{15625g}{720}=?`

we have put question mark on RHS of equation no 7 because we dont know whether their sum will be zero or no and morever we dont need that equation because we have six unknown so we want only 6 equations

converting it into matrix form for solving the unknowns

`[[1,1,1,1,1,1],[0,1,2,3,4,5],[0,1/2,2,9/2,8,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]]*[[a],[b],[c],[d],[e],[g]]=[[0],[0],[1],[0],[0],[0]]`

AX = B ........(6)

where A = coefficient matrix

X = unknown matrix

B = constant matrix

premultiplying by `A^-1`

`A^-1 AX = A^-1 B`

`IX=A^-1B`

`X=A^-1B`

multiplying the matrix B with inverse of matrix A

`[[a],[b],[c],[d],[e],[g]]=[[1,-2.28333,3.75,-4.25,3,-1],[0,5,-12.83333,17.75,-14,5],[0,-5,17.83333,-29.5,26,-10],[0,3.33333,-13,24.5,-24,10],[0,-1.25,5.08333,-10.25,11,-5],[0,0.2,-0.83333,1.75,-2,1]]*[[0],[0],[1],[0],[0],[0]]`

on solving the above matrix multiplication we get coefficients values as

a = 3.75 , b = -12.83333 , c = 17.83333 , d = -13 , e = 5.083333 , g = -0.83333

NOTE - For cross checking our coefficients are right, the sum of all the coefficient should be equal to zero

therefore our final formula for second derivative 4th order approximation is -

`frac{del^2u}{delx^2}=frac{3.75f(i)-12.83333f(i+1)+17.83333f(i+2)-13f(i+3)+5.08333f(i+4)-0.83333f(i+5)}{Deltax^2}`

now consider equation number 7

`0+b/720+frac{64c}{720}+frac{729d}{720}+frac{4096e}{720}+frac{15625g}{720}=0-12.83333/720+frac{64*17.83333}{720}-frac{729*13}{720}+frac{4096*5.08333}{720}-frac{15625*0.83333}{720}=-0.761058`

as we see that sum of equation number 7 is not zero hence in our final equation for approximation we will only have 2nd and 6th order derivative

***deriving the 4th order apprximation for skewed right hand side differencing***

`frac{del^2u}{delx^2} = (0+b/2+2c+frac{9d}{2}+frac{16e}{2}+frac{25g}{2})f^2(i) Deltax^2+ (0+b/720+frac{64c}{720}+frac{729d}{720}+frac{4096e}{720}+frac{15625g}{720})f^6(i) Deltax^6`

dividing the above equation by `Deltax^2`

`frac{del^2u}{delx^2} = (0+b/2+2c+frac{9d}{2}+frac{16e}{2}+frac{25g}{2})f^2(i)+ (0+b/720+frac{64c}{720}+frac{729d}{720}+frac{4096e}{720}+frac{15625g}{720})f^6(i) Deltax^4`

as we can see that the truncation error is of order 4,hence it is called as 4th order approximation

• SKEWED LEFT SIDE 4th ORDER APPROXIMATION

`frac{del^2u}{delx^2}=frac{af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)}{Deltax^2}`

drawing taylors table-

	`f(i)`	`Deltaxf'(i)`	`Deltax^2f''(i)`	`Deltax^3f'''(i)`	`Deltax^4f''''(i)`	`Deltax^5f''''(i)`	`Deltax^6f'''''(i)`
`af(i-5)`	a	-5a	25a/2	-125a/6	625a/24	-3125a/120	15625a/720
`bf(i-4)`	b	-4b	16b/2	-64b/6	256b/24	-1024b/120	4096b/720
`cf(i-3)`	c	-3c	9c/2	-27c/6	81c/24	-243c/120	729c/720
`df(i-2)`	d	-2d	4d/2	-8d/6	16d/24	-32d/120	64d/720
`ef(i-1)`	e	-e	e/2	-e/6	e/24	-e/120	e/720
`gf(i)`	g	0	0	0	0	0	0
sum of all the coefficients column vise `frac{del^2u}{delx^2}`	0	0	1	0	0	0	?

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

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

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

`-frac{125a}{6}-frac{64b}{6}-frac{27c}{6}-frac{8d}{6}-frac{e}{6}-0=0`

`frac{625a}{24}+frac{256b}{24}+frac{81c}{24}+frac{16d}{24}+e/24+0=0`

`-frac{3125a}{120}-frac{1024b}{120}-frac{243c}{120}-frac{32d}{120}-e/120-0=0`

`frac{15625a}{720}+frac{4096b}{720}+frac{729c}{720}+frac{64d}{720}+e/720+0=?`

we have put question mark on RHS of equation no 7 because we dont know whether their sum will be zero or no and morever we dont need that equation because we have six unknown so we want only 6 equations

converting it into matrix form for solving the unknowns

`[[1,1,1,1,1,1],[-5,-4,-3,-2,-1,0],[25/2,8,9/2,2,1/2,0],[-125/6,-64/6,-27/6,-8/6,-1/6,0],[625/24,256/24,81/24,16/24,1/24,0],[-3125/120,-1024/120,-243/120,-32/120,-1/120,0]]*[[a],[b],[c],[d],[e],[g]]=[[0],[0],[1],[0],[0],[0]]`

AX = B ........(6)

where A = coefficient matrix

X = unknown matrix

B = constant matrix

premultiplying by `A^-1`

`A^-1 AX = A^-1 B`

`IX=A^-1B`

`X=A^-1B`

multiplying the matrix B with inverse of matrix A

`[[a],[b],[c],[d],[e],[g]]=[[0,-0.2,-0.83333,-1.75,-2,-1],[0,1.25,5.08333,10.25,11,5],[0,-3.33333,-13,-24.5,-24,-10],[0,5,17.83333,29.5,26,10],[0,-5,-12.83333,-17.75,-24,-5],[1,2.28333,3.75,4.25,3,1]]*[[0],[0],[1],[0],[0],[0]]`

on solving the above matrix multiplication we get coefficients values as

a = -0.83333 , b = 5.08333 , c = -13 , d = 17.83333 , e = -12.83333 , g = 3.75

therefore our final formula for second derivative 4th order approximation is -

`frac{del^2u}{delx^2}=frac{-0.83333f(i-5)+5.08333(i-4)-13f(i-3)+17.83333f(i-2)-12.83333f(i-1)+3.75f(i)}{Deltax^2}`

now consider equation number 7

`frac{15625a}{720}+frac{4096b}{720}+frac{729c}{720}+frac{64d}{720}+frac{e}{720}+0=-frac{15625*0.83333}{720}+frac{4096*5.08333}{720}-frac{729*13}{720}+frac{64*17.83333}{720}-frac{12.83333}{720}+0=-0.761058`

as we see that sum of equation number 7 is not zero hence in our final equation for approximation we will only have 2nd and 6th order derivative

***deriving the 4th order apprximation for skewed left hand side differencing***

`frac{del^2u}{delx^2} = (frac{25a}{2}+8b+frac{9c}{2}+frac{4d}{2}+frac{e}{2}+0)f^2(i) Deltax^2+ (frac{15625a}{720}+frac{4096b}{720}+frac{729c}{720}+frac{64d}{720}+frac{e}{720})f^6(i) Deltax^6`

dividing the above equation by `Deltax^2`

`frac{del^2u}{delx^2} = (frac{25a}{2}+8b+frac{9c}{2}+frac{4d}{2}+frac{e}{2}+0)f^2(i)+ (frac{15625a}{720}+frac{4096b}{720}+frac{729c}{720}+frac{64d}{720}+frac{e}{720})f^6(i) Deltax^4`

as we can see that the truncation error is of order 4,hence it is called as 4th order approximation

The approximation is called 4th order because the order of the truncation error is 4

4. Code :

Main code - taylor_table_challenge

%solving a function f(x)=exp(x)* cos(x)
%solving the above function analytically and numerically
%f''(x) = -2*exp(x)*sin(x)

clear all
close all
clc


x=pi/2;

%solving for a range of dx
dx=linspace(pi/6,pi/6000,30);

for i = 1:length(dx)
  
  
  %calculating errors for different schemes
  error_central_differencing(i) = central_differencing(x,dx(i));
  
  error_skewed_right_differencing(i) = skewed_right_differencing(x,dx(i));
  
  error_skewed_left_differencing(i)  = skewed_left_differencing(x,dx(i));
  
  
end


%plotting of results
figure(1);
loglog(dx,error_central_differencing,'LineWidth',2,'color','r');
hold on
loglog(dx,error_skewed_right_differencing,'LineWidth',2,'color','g');
loglog(dx,error_skewed_left_differencing,'LineWidth',2,'color','k');
xlabel('Grid size (dx)');
ylabel('Error');
legend({'central differencing','skewed right differencing','skewed left differencing'},'Location','Northwest');

function calling code - central_differencing(x,dx)

function out = central_differencing(x,dx)
  
  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];
  
  B = [0;0;1;0;0];

  X = inv(A)*B;
  
  %solving analytical function exp(x)*cos(x)
  analytical = -2*exp(x)*sin(x);
  
  %central differencing of 4th order approximation
  central_differencing = ((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;
  
  %error
  out = abs(central_differencing - analytical);
  
  
endfunction

function calling code - skewed_right_differencing(x,dx)

function out = skewed_right_differencing(x,dx)
  
  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 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];
  
  B = [0;0;1;0;0;0];
  
  X = inv(A)*B;
  
  %solving analytical function exp(x)*cos(x)
  analytical = -2*exp(x)*sin(x);
  
  %skewed right side of 4th order approximation
  skewed_right_differencing = ((X(1)*exp(x)*cos(x)) + (X(2)*exp(x+dx)*cos(x+dx)) + (X(3)*exp(x+(2*dx))*cos(x+(2*dx))) + (X(4)*exp(x+(3*dx))*cos(x+(3*dx))) + (X(5)*exp(x+(4*dx))*cos(x+(4*dx))) + (X(6)*exp(x+(5*dx))*cos(x+(5*dx))))/dx^2;
  
  %error
  out = abs(skewed_right_differencing - analytical);
  
  end

function calling code - skewed_left_differencing(x,dx)

function out = skewed_left_differencing(x,dx)
  
  A = [1 1 1 1 1 1 ; -5 -4 -3 -2 -1 0 ; 25/2 8 9/2 2 1/2 0 ; -125/6 -64/6 -27/6 -8/6 -1/6 0 ; 625/24 256/24 81/24 16/24 1/24 0 ; -3125/120 -1024/120 -243/120 -32/120 -1/120 0];
  
  B = [0;0;1;0;0;0];
  
  X = inv(A)*B;
  
  %solving analytical function exp(x)*cos(x)
  analytical = -2*exp(x)*sin(x);
  
  %skewed left side of 4th order approximation
  skewed_left_differencing = ((X(1)*exp(x-5*dx)*cos(x-5*dx)) + (X(2)*exp(x-4*dx)*cos(x-4*dx)) + (X(3)*exp(x-3*dx)*cos(x-3*dx)) + (X(4)*exp(x-2*dx)*cos(x-2*dx)) + (X(5)*exp(x-dx)*cos(x-dx)) + (X(6)*exp(x)*cos(x)))/dx^2;
  
  %error
  out = abs(skewed_left_differencing - analytical);
  
  
  end

5. Results :

6. Conclusion : The following points we can conclude from the above results

• The error for central differencing scheme is less than both skewed left and right side differencing scheme
• The error for skewed right and left side differencing is almost same
• The error value decreases as we decrease the grid spacing (dx)

Advantages of Skewed differencing scheme over Central differencing

• For central differencing scheme we should have information of the neighbouring points i.e both left and right points whereas for skewed differencing scheme we want either information of left side points for skewed left differencing scheme or information of right side points for skewed right differencing scheme
• For heat transfer problem in 1D plate flowing from left to right, we cannot use central differencing scheme at boundary because we dont have any points left to it but we can use central differencing scheme in middle of the plate hence we can use skewed differencing scheme at the boundary and central differencing scheme in middle of the discrete points for more accurate results.