Week 3 - Taylor table method and Matlab code

Our aim for the challenge is to derive 4th order approximation of 2nd order derivative using central difference, skewed right sided difference and skewed left sided difference. we also need to find 2nd order derivative of analytical function exp(x)*cos(x).

Discritization:

Numerical discritization is the technique which is used to convert an PDE into algebric equation which can be computed using programming language. This technique employes concepts of grid points and grid spacing using taylors tabel expansion.

Taylor's expansion:

`f(x+dx)=f(x)+f'(x)*dx/1!+f''(x)*(dx)^2/2!+f'''(x)*(dx)^3/3!+......+f^n(x)*(dx)^n/n!`

where dx is the grid spacing over a numerical stencil as below

The method is employed in order to get more accurate and refined mesh higher the number of approximation lower the mesh point and it increases the computational speed.

1)Central difference:

let us consider the central difference scheme for which we wil need a grid point at the centre f(i) and grid points on a particular grid spacing dx on both side of that point.

lets us consider 5 grid points with x as the centre.

f(x-2dx),f(x-dx),f(x),f(x+dx),f(x+2dx)

therefore we approximate the 2nd order derivative as

`(d^2f)/dx^2=a*f(x-2dx)+b*(x-dx)+c*(x)+d*f(x+dx)+e*f(x+2dx)`

To find the 4th order approximation of 2nd order derivative we have to take coefficient of 2nd order 1.

the taylor tabel is as follows:

The tabel gives us 5 eqations and 5 unknowns which we can solve using matrix inversion in matlab.

a+b+c+d+e=0

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

a+b+0+d+4e=1

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

16a+b+0+d+16e=0

on solving the equation we get a=-0.0833 b=1.3333 c=-2.5 d=1.3333 e=-0.0833

as a=e and b=d the f^5(x)*(Δ x)^5/5! = 0

and f^6(x)*(Δ x)^6/6! is approximation term hence we get 4th order approximation of 2nd order diffrentiation as:

(d^2f)/dx^2=-0.08333f(x-2dx)+1.333*(x-dx)-2.5*(x)+1.333*f(x+dx)-0.08333*f(x+2dx)/dx^2

2) 4th order skewed right sided 2nd order approximation:

Let us consider 5 grid points for skewed right sided scheme

f(x),f(x+dx),f(x+2dx),f(x+dx),f(x+4dx)

approximation for 2nd order derivative of skewed right sided :

(d^2f)/dx^2=a*f(x)+b*(x+dx)+c*(x+2dx)+d*f(x+3dx)+e*f(x+4dx)

To find the 4th order approximation of 2nd order derivative we have to take coefficient of 2nd order 1.

the taylor tabel is as follows:

The tabel gives us 5 eqations and 5 unknowns which we can solve using matrix inversion in matlab.

a+b+c+d+e=0

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

0+b+4c+9d+16e=1

0+b+8c+27d+64e=0

0+b+16c+81d+256e=0

on solving the equation we get a=2.9167 b=-8.6667 c=9.5 d=-4.6667 e=0.91667

as approximately a=e and b=d the f^5(x)*(Δ x)^5/5! = 0

and f^6(x)*(Δ x)^6/6! is approximation term hence we get 4th order approximation of 2nd order diffrentiation as:

(d^2f)/dx^2=2.9167*f(x)-8.6667*(x+dx)+9.5*(x+2dx)-4.6667*f(x+3dx)+0.91667*f(x+dx)/dx^2.

3) 4th order skewed right sided 2nd order approximation:

Let us consider 5 grid points for skewed right sided scheme

f(x),f(x+dx),f(x+2dx),f(x+dx),f(x+4dx)

approximation for 2nd order derivative of skewed right sided :

(d^2f)/dx^2=a*f(x)+b*(x-dx)+c*(x-2dx)+d*f(x-3dx)+e*f(x-4dx)

To find the 4th order approximation of 2nd order derivative we have to take coefficient of 2nd order 1.

the taylor tabel is as follows:

The tabel gives us 5 eqations and 5 unknowns which we can solve using matrix inversion in matlab.

a+b+c+d+e=0

0-b+-2c-3d-4e=0

0+b+4c+9d+16e=1

0-b-8c-27d-64e=0

0+b+16c+81d+256e=0

on solving the equation we get a=2.9167 b=-8.6667 c=9.5 d=-4.6667 e=0.91667

as approximately a=e and b=d the f^5(x)*(Δ x)^5/5! = 0

and f^6(x)*(Δ x)^6/6! is approximation term hence we get 4th order approximation of 2nd order diffrentiation as:

(d^2f)/dx^2=2.9167*f(x)-8.6667*(x+dx)+9.5*(x+2dx)-4.6667*f(x+3dx)+0.91667*f(x+dx)/dx^2.

Matlab Code:

clear all
close all
clc
%input variable
x=pi/3;
%range of dx values
dx=linspace(pi/3,pi/400,10);
%for loop to compute value of 2nf order diffrentiation with 4th order approximation 
for i=1:length(dx)
  central_error(i)=central_difference(x,dx(i));
  skewed_right_error(i)=skewed_right(x,dx(i));
  skewed_left_error(i)=skewed_left(x,dx(i));
endfor
%plot 
figure(1)
plot(dx,central_error,'color','b','marker','*');
hold on
plot(dx,skewed_right_error,'color','r','marker','+');
hold on 
plot(dx,skewed_left_error,'color','y','marker','d');
grid on 
legend('centralerror','skewedrighterror','skewedlefterror')
xlabel('dx')
ylabel('errors')

The above code computes the value of error for range of dx using a for loop. it stores the value of each central error, skewed right error and skewed left error by calling function for central_diffrence,skewed_right and skewed left.

The function computes the abs error value between analytical value and approximation values

the analytical second order derivative for function is -2*exp(x)*sin(x)

Function code:

1)Central Difference:

function output=central_difference(x,dx)
  analytical=-2*exp(x)*sin(x);
  central_difference=(-1*exp(x-2*dx)*cos(x-2*dx)+16*exp(x-dx)*cos(x-dx)-30*exp(x)*cos(x)+16*exp(x+dx)*cos(x+dx)-1*exp(x+2*dx)*cos(x+2*dx))/12*dx^2
output=abs(analytical-central_difference)
  endfunction

2)Skewed right-sided difference:

function output=skewed_right(x,dx)
  analytical=-2*exp(x)*sin(x);
  skewed_right_sided_diff=((2.9167*exp(x)*cos(x))+(-8.6667*exp(x+dx)*cos(x+dx))+(9.5*exp(x+2*dx)*cos(x+2*dx))+(-4.6667*exp(x+3*dx)*cos(x+3*dx))+(0.91667*exp(x+4*dx)*cos(x+4*dx)))/dx^2
  output=abs(analytical-skewed_right_sided_diff)
endfunction

3)Skewed left-sided difference:

function output=skewed_left(x,dx)
  analytical=-2*exp(x)*sin(x);
  skewed_left_sided_diff=((2.9167*exp(x)*cos(x))+(-8.6667*exp(x-dx)*cos(x-dx))+(9.5*exp(x-2*dx)*cos(x-2*dx))+(-4.6667*exp(x-3*dx)*cos(x-3*dx))+(0.91667*exp(x-4*dx)*cos(x-4*dx)))/dx^2
  output=abs(analytical-skewed_left_sided_diff)
endfunction

PLOT:

Conclusion:

The plot showns that the error is very less for central difference method as comapred to skewed right sided and skewed left sided.the error in central difference is less as we consider information on both the sides of numerical stencil.

But the reason we need skewed right sided and skewed left sided method is that when we need to apply approximation of boundary conditions we do not have grid points on both the side hence in such conditions we need to use skewed methods.