Week 3 - Taylor table method and Matlab code

TYLOR TABLE -Method and MATLAB code                Date: - 19-01-2021

 

4^th order approximations of the second-order derivative: -

• Our main aim to find the difference or absolute error between exact value and approximate values from three approximations.
• Make a plot between ‘dx’ and ‘three absolute errors’.

 

Exact value: -

Given function                                    f(x) = y = e^(x)* cos(x)

                                    dy/dx = e^(x)(-sin(x)) + (e^(x)*cos(x))

                                    d²y/dx² = -2*sin(x)*e^(x)

• Consider x =л/3 and dx =50 values between 0.1 to 0.01

d²y/dx²|_x=л/3 = -4.935745353035276 = Exact Value

Approximate values: -

• We have to find three 4^th order approximations for second order derivative of above given function.
• Derivations of the central difference skewed right side difference and skewed left side difference approximations are mentioned below.

4^th order approximation of Central difference: -

• For central difference scheme,

No of Nodes =  p+q-1

                        Where,

                        p – order of derivative

                        q – order of approximation

• So, No of nodes = 2+4-1 = 5

Δx²(ꝺ²f/ꝺx²) = a*f(i-2) + b*f(i-1) + c*f(i) + d*f(i+1) + e*f(i+2)

Where, a,b,c,d,e are coefficients

• Here our main work is to find coeffects by using Tylor’s table

 

 

 

Δx2(ꝺ2f/ꝺx2)	f(i)	Δx*f'(i)	Δx^2*f''(i)	Δx^3*f'''(i)	Δx^4*f''''(i)
a*f(i-2)	a	-2a	2a	-4a/3	2a/3
b*f(i-1)	b	-b	b/2	-b/6	b/24
c*f(i)	c	0	0	0	0
d*f(i+1)	d	d	d/2	d/6	d/24
e*f(i+2)	e	2e	4e/2	4e/3	2e/3
sum	0	0	1	0	0

Table 1:- 2^nd order derivative Tylor’s table for given function(central difference)

• Now we can create MATLAB code for finding coefficients.
• Creating matrix in MATLAB as using following code

Code:-

%% central diffrence for 4th order approximation

u1=zeros(5);

u1(1,:)=[1,1,1,1,1];

u1(2,:)=[-2,-1,0,1,2];

u1(3,:)=[2,1/2,0,1/2,2];

u1(4,:)=[-4/3,-1/6,0,1/6,4/3];

u1(5,:)=[2/3,1/24,0,1/24,2/3];

v1=zeros(5,1);

v1(3)=1;

w1=inv(u1)*v1;

 

 

• From above code we can have coefficients values as follows

 

W₁  =

    

a	-0.08333
b	1.333333
c	-2.5
d	1.333333
e	-0.08333

 

 

 

 

 

 

4^th order approximation of Skewed right sided difference: -

 

• For Skewed right sided scheme,

No of Nodes =  p+q

                        Where,

                        p – order of derivative

                        q – order of approximation

• So, No of nodes = 2+4 = 6

Δx²(ꝺ²f/ꝺx²) = a*f(i) + b*f(i+1) + c*f(i+2) + d*f(i+3) + e*f(i+4) + g*f(i+5) ------------- 1

Where, a,b,c,d,e,g are coefficients

• Here our main work is to find coeffects by using Tylor’s table

 

Δx2(ꝺ2f/ꝺx2)	f(i)	Δx*f'(i)	Δx^2*f''(i)	Δx^3*f'''(i)	Δx^4*f''''(i)	Δx^5*f'''''(i)
a*f(i)	a	0	0	0	0	0
b*f(i+1)	b	b	b/2	b/6	b/24	b/120
c*f(i+2)	c	2c	2c	4c/3	2c/3	32b/120
d*f(i+3)	d	3d	9d/2	27d/6	81d/24	243d/120
e*f(i+4)	e	4e	16e/2	64e/6	256e/24	1024e/120
g*f(i+5)	g	5g	25g/2	125g/6	625g/24	3125g/120
sum	0	0	1	0	0	0

Table 2: - 2^nd order derivative Tylor’s table for given function (skewed right side)

 

• Now we can create MATLAB code for finding coefficients.
• Creating matrix in MATLAB as using following code

 

Code:

%% Skewed RSD

u2=zeros(6);

u2(1,:)=[1,1,1,1,1,1];

u2(2,:)=[0,1,2,3,4,5];

u2(3,:)=[0,1/2,2,9/2,8,25/2];

u2(4,:)=[0,1/6,4/3,27/6,64/6,125/6];

u2(5,:)=[0,1/24,2/3,81/24,256/24,625/24];

u2(6,:)=[0,1/120,32/120,243/120,1024/120,3125/120];

v2=zeros(6,1);

v2(3)=1;

w2=inv(u2)*v2;

 

• From above code we can have coefficients values as follows

 

W₂  =

a	3.75
b	-12.8333
c	17.83333
d	-13
e	5.083333
g	-0.83333

 

4^th order approximation of Skewed left sided difference: -

 

• For Skewed right sided scheme,

No of Nodes =  p+q

                        Where,

                        p – order of derivative

                        q – order of approximation

• So, No of nodes = 2+4 = 6

Δx²(ꝺ²f/ꝺx²) = a*f(i-5) + b*f(i-4) + c*f(i-3) + d*f(i-2) + e*f(i-1) + g*f(i)  ---------- 2

Where, a,b,c,d,e,g are coefficients

• Here our main work is to find coeffects by using Tylor’s table

 

Δx2(ꝺ2f/ꝺx2)	f(i)	Δx*f'(i)	Δx^2*f''(i)	Δx^3*f'''(i)	Δx^4*f''''(i)	Δx^5*f'''''(i)
a*f(i-5)	a	-5a	25a/2	-125a/6	625a/24	-3125a/120
b*f(i-4)	b	-4b	16b/2	-64b/6	256b/24	-1024b/120
c*f(i-3)	c	-3c	9c/2	-27c/6	81c/24	-243c/120
d*f(i-2)	d	-2d	2d	-4d/3	2d/3	-32d/20
e*f(i-1)	e	-e	e/2	-e/6	e/24	-e/120
g*f(i)	g	0	0	0	0	0
sum	0	0	1	0	0	0

Table 2: - 2^nd order derivative Tylor’s table for given function (skewed left side)

 

 

 

• Now we can create MATLAB code for finding coefficients.
• Creating matrix in MATLAB as using following code

 

Code:

%% Skewed LSD

u3=zeros(6);

u3(1,:)=[1,1,1,1,1,1];

u3(2,:)=[-5,-4,-3,-2,-1,0];

u3(3,:)=[25/2,8,9/2,2,1/2,0];

u3(4,:)=[-125/6,-64/6,-27/6,-4/3,-1/6,0];

u3(5,:)=[625/24,256/24,81/24,2/3,1/24,0];

u3(6,:)=[-3125/120,-1024/120,-243/120,-32/120,-1/120,0];

v3=zeros(6,1);

v3(3)=1;

w3=inv(u3)*v3;

 

• From above code we can have coefficients values as follows

                           

W₃  =

a	-0.83333
b	5.083333
c	-13
d	17.83333
e	-12.8333
g	3.75

 

Steps to find desire output: -

• By using above coefficients, we can find output by using for loop
• Here we are using calling function for given ‘y’ function

Function code: -

function y=numerical(x)

 

        y =exp(x)*cos(x);

 

end

Code: -

%%

x=pi/3;

dx=linspace(0.1,0.01,50);

exact_value = -2*sin(x)*exp(x);

for i=1:length(dx)

f_central(i)=[w1(1)*numerical(x-2*dx(i))+w1(2)*numerical(x-dx(i))+w1(3)*numerical(x)+w1(4)*numerical(x+dx(i))+w1(5)*numerical(x+2*dx(i))]/(dx(i)^2);

f_RSD(i)=[w2(1)*numerical(x)+w2(2)*numerical(x+dx(i))+w2(3)*numerical(x+2*dx(i))+w2(4)*numerical(x+3*dx(i))+w2(5)*numerical(x+4*dx(i))+w2(6)*numerical(x+5*dx(i))]/(dx(i)^2);

f_LSD(i)=[w3(1)*numerical(x-5*dx(i))+w3(2)*numerical(x-4*dx(i))+w3(3)*numerical(x-3*dx(i))+w3(4)*numerical(x-2*dx(i))+w3(5)*numerical(x-dx(i))+w3(6)*numerical(x)]/(dx(i)^2);

a(i)=abs(f_central(i)-exact_value);

b(i)=abs(f_RSD(i)-exact_value);

c(i)=abs(f_LSD(i)-exact_value);

end

 

• This output truncation error is very small so, we are using log graph to make plot visible.

Code:

figure(1)

loglog(dx,a,'-r','linewidth',3);

hold on;

loglog(dx,b,'-k','linewidth',3);

hold on;

loglog(dx,c,'-b','linewidth',3);

grid on;

xlabel('dx');

ylabel('Absolute error');

legend('Central Difference of 4th Order','Skewed Right Side Difference of 4th Order','Skewed Left Side Difference of 4th Order');

title('Consistency of different scheme');

 

OUTPUT GRAPH: -

 

 

Proof of Skewed schemes are 4^th order approximation: -

Skewed right sided difference: -

• From equation 1 we can write that

 

Δx²(ꝺ²f/ꝺx²) = a*f(i) + b*f(i+1) + c*f(i+2) + d*f(i+3) + e*f(i+4) + g*f(i+5)

• Then expand functions by using Tylor series.

L.H.S =    a*f(i) + b*[f(i) + Δx*f’(i) + Δx² * f’’(i)/2! + Δx³ * f’’’(i)/3! + Δx⁴ * f’’’’(i)/4! + Δx⁵ *  f’’’’’(i) + Δx⁶ * f’’’’’’(i)/6! ………..etc] + c*[f(i) + 2Δx*f’(i) + (2Δx)² * f’’(i)/2! + (2Δx)³ * f’’’(i)/3! + (2Δx)⁴ * f’’’’(i)/4! + (2Δx)⁵ * f’’’’’(i) + (2Δx)⁶ * f’’’’’’(i)/6! ……… etc] + d*[f(i) + 3Δx*f’(i) + (3Δx)² * f’’(i)/2! + (3Δx)³ * f’’’(i)/3! + (3Δx)⁴ * f’’’’(i)/4! + (3Δx)⁵ * f’’’’’(i) + (3Δx)⁶ * f’’’’’’(i)/6! ……… etc] + e*[f(i) + 4Δx*f’(i) + (4Δx)² * f’’(i)/2! + (4Δx)³ * f’’’(i)/3! + (4Δx)⁴ * f’’’’(i)/4! + (4Δx)⁵ * f’’’’’(i) + (4Δx)⁶ * f’’’’’’(i)/6! ……… etc] + g*[f(i) + 5Δx*f’(i) + (5Δx)² * f’’(i)/2! + (5Δx)³ * f’’’(i)/3! + (5Δx)⁴ * f’’’’(i)/4! + (5Δx)⁵ * f’’’’’(i) + (5Δx)⁶ * f’’’’’’(i)/6! ……… etc]

LHS = (a+b+c+d+e+g) * f(i) + (b+2c+3d+4e+5g) *Δx *f’(i) + (b/2+2c+9d/2+8e+25g/2) *Δx² *f’’(i) + (b/6+8c/6+27d/6+64e/6+125g/6) *Δx³ *f’’’(i) + (b+16c+81d+256e+625g) *(Δx⁴/24)*f’’’’(i) + (b+32c+243d+1024e+3125g) *(Δx⁵/120) *f’’’’’(i) + (b+64c+729d+4096e+15625g) *(Δx⁶/720) *f’’’’’’(i) + ……………..  etc

• Substitute values of coefficients a=3.7 , b=-12.8333 , c=17.8333 , d=-13 , e= 5/08333 and g= -0.8333 then we endup with

LHS = Δx²(ꝺ²f/ꝺx²) =  I * Δx² *f’’(i) + J * Δx⁶ * f’’’’’’(i)

Where I and J are the constant non zero values after substituting the values of coefficients .

 

ꝺ²f/ꝺx² = =  I *f’’(i) + J * Δx⁴ * f’’’’’’(i)  

 

• From above equation we can say that order of the 2^nd order derivative for skewed right sided difference is 4^th
• Hence proved ………..

 

Skewed Left sided difference: -

• From equation 2 we can write that

 

Δx²(ꝺ²f/ꝺx²) = a*f(i-5) + b*f(i-4) + c*f(i-3) + d*f(i-2) + e*f(i-1) + g*f(i)

• Then expand functions by using Tylor series.

L.H.S =    g*f(i) + e*[f(i) -  Δx*f’(i) + Δx² * f’’(i)/2! - Δx³ * f’’’(i)/3! + Δx⁴ * f’’’’(i)/4! - Δx⁵ *  f’’’’’(i) + Δx⁶ * f’’’’’’(i)/6! ………..etc] + d*[f(i) - 2Δx*f’(i) + (2Δx)² * f’’(i)/2! - (2Δx)³ * f’’’(i)/3! + (2Δx)⁴ * f’’’’(i)/4! - (2Δx)⁵ * f’’’’’(i) + (2Δx)⁶ * f’’’’’’(i)/6! ……… etc] + c*[f(i) - 3Δx*f’(i) + (3Δx)² * f’’(i)/2! - (3Δx)³ * f’’’(i)/3! + (3Δx)⁴ * f’’’’(i)/4! - (3Δx)⁵ * f’’’’’(i) + (3Δx)⁶ * f’’’’’’(i)/6! ……… etc] + b*[f(i) - 4Δx*f’(i) + (4Δx)² * f’’(i)/2! - (4Δx)³ * f’’’(i)/3! + (4Δx)⁴ * f’’’’(i)/4! - (4Δx)⁵ * f’’’’’(i) + (4Δx)⁶ * f’’’’’’(i)/6! ……… etc] + a*[f(i) - 5Δx*f’(i) + (5Δx)² * f’’(i)/2! - (5Δx)³ * f’’’(i)/3! + (5Δx)⁴ * f’’’’(i)/4! - (5Δx)⁵ * f’’’’’(i) + (5Δx)⁶ * f’’’’’’(i)/6! ……… etc]

• Similary as skewed right sided difference , Substitute values of coefficients g=3.7 , e=-12.8333 , d=17.8333 , c=-13 , b= 5/08333 and a= -0.8333 then we endup with

LHS = Δx²(ꝺ²f/ꝺx²) =  I * Δx² *f’’(i) + J * Δx⁶ * f’’’’’’(i)

Where I and J are the constant non zero values after substituting the values of coefficients .

 

ꝺ²f/ꝺx² = =  I *f’’(i) + J * Δx⁴ * f’’’’’’(i)  

 

• From above equation we can say that order of the 2^nd order derivative for skewed left sided difference is 4^th
• Hence proved ………..

 

 

Advantages of Skewed schemes over Central difference schemes: -

• In central difference we must have information of symmetric points right side and left side of desire point (point at which we need information). But in skewed schemes it is not important to be symmetric.

These skewed schemes are very important at finding the first point information means after boundary condition point. In figure consider domain and ‘A’, ‘B’ are boundary point where we have information. Then, by using skewed schemes we can able to find information at point ‘B’ and ‘G’ very easily but with central difference we can’t.

%%Function programme
function y=numerical(x)

        y =exp(x)*cos(x);%given function in problem

end

%%function file and main programme file should be in same folder

%%main programme
clear all;
close all;
clc;
%% central diffrence for 4th order approximation
u1=zeros(5);%defining null matrix of size 5X5
u1(1,:)=[1,1,1,1,1];
u1(2,:)=[-2,-1,0,1,2];
u1(3,:)=[2,1/2,0,1/2,2];
u1(4,:)=[-4/3,-1/6,0,1/6,4/3];
u1(5,:)=[2/3,1/24,0,1/24,2/3];
v1=zeros(5,1);
v1(3)=1;
w1=inv(u1)*v1;%by using inverse operator we can able to find coefficients 
%% Skewed RSD
u2=zeros(6);%defining null matrix of size 6X6
u2(1,:)=[1,1,1,1,1,1];
u2(2,:)=[0,1,2,3,4,5];
u2(3,:)=[0,1/2,2,9/2,8,25/2];
u2(4,:)=[0,1/6,4/3,27/6,64/6,125/6];
u2(5,:)=[0,1/24,2/3,81/24,256/24,625/24];
u2(6,:)=[0,1/120,32/120,243/120,1024/120,3125/120];
v2=zeros(6,1);
v2(3)=1;
w2=inv(u2)*v2;%by using inverse operator we can able to find coefficients
 
%%
x=pi/3;
dx=linspace(0.1,0.01,50);
exact_value = -2*sin(x)*exp(x);
for i=1:length(dx)
f_central(i)=[w1(1)*numerical(x-2*dx(i))+w1(2)*numerical(x-dx(i))+w1(3)*numerical(x)+w1(4)*numerical(x+dx(i))+w1(5)*numerical(x+2*dx(i))]/(dx(i)^2);%assigning values of coefficients and finding truncation error
f_RSD(i)=[w2(1)*numerical(x)+w2(2)*numerical(x+dx(i))+w2(3)*numerical(x+2*dx(i))+w2(4)*numerical(x+3*dx(i))+w2(5)*numerical(x+4*dx(i))+w2(6)*numerical(x+5*dx(i))]/(dx(i)^2);%assigning values of coefficients and finding truncation error
f_LSD(i)=[w3(1)*numerical(x-5*dx(i))+w3(2)*numerical(x-4*dx(i))+w3(3)*numerical(x-3*dx(i))+w3(4)*numerical(x-2*dx(i))+w3(5)*numerical(x-dx(i))+w3(6)*numerical(x)]/(dx(i)^2);%assigning values of coefficients and finding truncation error
a(i)=abs(f_central(i)-exact_value);
b(i)=abs(f_RSD(i)-exact_value);
c(i)=abs(f_LSD(i)-exact_value);
end
%%

figure(1)
loglog(dx,a,'-r','linewidth',3);%%using log plot so that we can able to see very minimum differenece in residuals
hold on;
loglog(dx,b,'-k','linewidth',3);
hold on;
loglog(dx,c,'-b','linewidth',3);
grid on;
xlabel('dx');
ylabel('Absolute error');
legend('Central Difference of 4th Order','Skewed Right Side Difference of 4th Order','Skewed Left Side Difference of 4th Order');
title('Consistency of different scheme');
%%