Week 3 - Taylor table method and Matlab code

Objective :- To derive the 4th order approximations of the second order derivative using

1. Central Difference

2. Skewed Right Sided Differene

3. Skewed Left Sided Difference

Introduction :- Taylor Table method is used to derive the equations for second order derivative using fourth order approximation by different schemes .

For Central Difference Scheme

No of nodes required = p+q-1 , where p is the order of the derivative and q is the order of approximation .So 2nd order derivative using 4th order approximation will require 2+4-1=5 no of nodes . So here 5 stencils are used for approximation .

Let `(d^2 f)/(dx^2)` be the equation to be numerically discretized .

Based on the above diagram the discretized equation is as follows :-

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

Using the taylor table method we determine the coefficients .

Taylor Table :-

We have five unknowns ,so we will use five equations to solve the problem . So the equations are :-

a+b+c+d+e=0

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

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

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

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

Solving the above equations by using inverse of the coefficient matrix in matlab we get coefficents as:-

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

Substituting the above values we have 6th order terms in R.H.S and 2nd order terms in L.H.S ,so dividing by 2nd order in R.H.S we get approximation of order 4 or we can say truncation error is of order 4 .

For Skewed Right Sided Difference Scheme

No of nodes required =p+q=2+4=6 .So here 6 no of nodes are required which is also the number of nodes required if you are calculating the same order derivative and same order approximation from left hand side. It is called skewed because there is no symmetry i.e the same no of nodes are not present on both left and right hand side of the current node .In the skewed right hand side the stencils take are on the right hand side of the current node and no stencil is present on the left side .

Based on the above diagram the discretized equation is as follows :-

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

Again using the taylor table method we get :-

We have 6 unkowns ,so 6 equations are required to solve the problem .So the equations are :-

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

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

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

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

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

b/120+4c/15+81d/40+128e/15+625g/24=0

Solving the above equations by same procedure in matlab we get :-

a=3.75 ,b=-12.83,c=17.83,d=-13,e=5.083,g=-0.8333

Substituting the values we get :-

`(d^2f)/dx^2=(3.75f(i)-12.83f(i+1)+17.83f(i+2)-13f(i+3)+5.083f(i+4)-0.8333f(i+5))/(trianglex^2)`

For Skewed Left Sided Difference Scheme

No of stencils =6

The discretized equation will be as follows :-

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

Using the taylor table method we get :-

Here also 6 unknowns , so 6 equations are required to solve the coefficients . The equations are :-

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

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

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

125a/6+32b/3+9c/2+4d/3+e/6=0

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

625a/24+128b/15+81c/40+4d/15+e/120=0

Solving the equations we get the coeffients as:-

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

Substituting the values we get

`(d^2f)/dx^2=(-0.8333f(i-5)+5.0833f(i-4)-13.0000f(i-3)+17.8333f(i-2)-12.8333f(i-1)+3.75f(i))/(trianglex^2)`

Matlab code for the evaluation of second order derivative by using different schemes :-

clc 
clear all
close all
x=50;
dx=linspace(10^(-5),10^(-1),100);
analytical_derivative =-2*(exp(x)*(sin(x)));
for i=1:length(dx)
   central_difference_approx(i)=((-5/60)*exp(x-2*dx(i))*cos(x-2*dx(i))+(4/3)*exp(x-dx(i))*cos(x-dx(i))-(10/4)*exp(x)*cos(x)+(4/3)*exp(x+dx(i))*cos(x+dx(i))+(-5/60)*exp(x+2*dx(i))*cos(x+2*dx(i)))/(dx(i)^2);
   central_difference_error(i)=abs(central_difference_approx(i)-analytical_derivative);
   skewed_right_hand_approx(i)=((15/4)*exp(x)*cos(x)-(77/6)*exp(x+dx(i))*cos(x+dx(i))+(107/6)*exp(x+2*dx(i))*cos(x+2*dx(i))-(13/1)*exp(x+3*dx(i))*cos(x+3*dx(i))+(61/12)*exp(x+4*dx(i))*cos(x+4*dx(i))-(5/6)*exp(x+5*dx(i))*cos(x+5*dx(i)))/(dx(i)^2);
   right_hand_error(i)=abs(skewed_right_hand_approx(i)-analytical_derivative);
   skewed_left_hand_approx(i)=((-5/6)*exp(x-5*dx(i))*cos(x-5*dx(i))+(61/12)*exp(x-4*dx(i))*cos(x-4*dx(i))-(13/1)*exp(x-3*dx(i))*cos(x-3*dx(i))+(107/6)*exp(x-2*dx(i))*cos(x-2*dx(i))+(-77/6)*exp(x-dx(i))*cos(x-dx(i))+(15/4)*exp(x)*cos(x))/(dx(i)^2);
   left_hand_error(i)=abs(skewed_left_hand_approx(i)-analytical_derivative);    
end
loglog(dx,central_difference_error,'linewidth',3)
hold on 
loglog(dx,right_hand_error,'linewidth',3)
hold on 
loglog(dx,left_hand_error,'linewidth',3)
xlabel('dx')
ylabel('Absolute error')
legend('Central Differencing Error ','Right Hand Skewed Error','Left Hand Skewed Error')
title('Error analysis of different schemes')

Output :-

For 100 nodes

For 500 nodes

a. Comparison of absolute error of different schemes :-

From the above plot ,it is found that right hand skewed rule gives approximately the same value as left hand skewed rule .But central differencing scheme is more accurate and preferable and the absolute error is less due to using information equally from both left hand side and from right hand side .Central differencing scheme otherwise uses this symmetry property from both nodes to evalaute the property at the current node .

b. Why a skewed scheme is useful ? What can a skewed scheme do that a CD scheme cannot do ?

Ans :- A skewed sytem may have higher truncation error in comparison to central differencing scheme . Both right and left skewed system have nearly same truncation error .But a skewed scheme is useful at all locations since it uses the knowlege of one node to evalaute the property at the current node .This is particularly helpful in evaluating the properties at the boundary .

Central Differencing Scheme whereas requires the knowledge of both left and right node .So this scheme is not useful in evaluating properties at the boundary or where a single node or stencil is present .If equivalent number of left and right nodes are not present then using central differencing scheme will be difficult since large error will creep up . Whereas skewed method can be used at any location if there is one node available for calculation .

Code Explaination :-

Step 1:- Clear the command window and work space by using commands clear all and close all .

Step 2:- The given function is exp(x)*cos(x). Find its derivative analytically . To find the second order derivative at a particular point we have to initiliaze so x is taken as 50 .x can be any arbitrary value .Numerical methods are used to evaluate so dx (small increment is taken in the range of (10^-5) to (10^-1) with 100 nodes and 500 nodes to study the behaviour of different schemes .

Step 3:- For loop is used here to calculate the second order derivative using central difference scheme ,skewed right hand and skewed left hand approximation for each node .Before that all coefficients are determined by matrix inversion method in matlab . Careful has to be taken in converting decimal into fraction by using sym command in matlab otherwise large truncation error will come into picture .

Step 4:- For plotting the graph use loglog to convert the values in 10powers so that the graph is meaningful .Hold on command is used so that in the same grapg other two methods variation can also be plotted .Linewidth is provided so that the lines are thick and the graph is visible .

Step 5:- x axis ,y axis label are provided using label command and what each line represent is given by legend command and title of the graph should explain the main theme of the graph using the command title .