Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method

Aim: To derive and write a program for in MATLAB for 4th order approximation for second order derivative by following methods:

1. Central difference

2. Skewed right sided difference

3. Skewed left sided difference

Objective:

1. To evaluate the second order derivative of the function f(x)=exp(x)*cosx, using above mentioned order of approximations.

2. To plot and compare the values and errors of the 3 methods derived.

Theory and derivations:

During solving a CFD Problem we may come across partial derivative equation which consists 2nd order derivative. Since due to complexity of the equation we may have to chose numerical approach rather than analytical. But when dealing with numerical approach we try to make the error margin as low as possible.

As the order of approximation get higher, we get higher accuracy. Hence, we are going to derive the three methods of approximation for an analytical function. The methods are following.

1. Central difference:

Here the nodes selected are symmetrical i.e the data is taken from both sides of point of focus. Using taylor's method, we get

`(del^2f)/(delx^2) = af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)`

Here there are two nodes taken from both side of the point of focus i.e f(i).

(NOTE: To select the number of nodes in consideration, we use following expression:

N=p+q-1 (For Central Differencing)

N=p+q (For skewed differencing)

Here N is no. of nodes, p is order of derivative and q is order of approximation.)

Taylor's Table for Central difference:

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

`[[1,1,1,1,1],[-2,-1,0,1,2],[4/2,1/2,0,1/2,4/2],[1,2,3,4,5],[-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]]`

on further solving we get the values of the unknown

i.e values a=-0.0833 b=1.333 c=-2.5000 d=1.333 e=-0.0833

These coefficients will further be used in code to get the value of the derivatibe.

1. Skewed right sided difference :- When the data point of the left nodes are not available ( at initital boundary conditions) then we have to rely on the available data set.

`(del^2f)/(delx^2) = af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+ff(i+5)`

Taylor's Table for Skewed right sided difference

	f(x)	dxf'(x)	dx^2f''(x)	dx^3f'''(x)	dx^4f''''(x)	dx^5f'''''(x)
af(i)	a	0	0	0	0	0
bf(i+1)	b	b	b/2	b/6	b/24	b/120
cf(i+2)	c	2c	4c/2	8c/6	16c/24	32c/120
df(i+3)	d	3d	9d/2	27d/6	81d/24	243d/120
ef(i+4)	e	4e	16e/2	64e/6	256e/24	1024e/120
ff(i+5)	f	5f	25f/2	125f/6	625f/24	3125f/120
af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+ff(i+5)	0	0	1	0	0	0

This table transform into matrix of type AX=B

`[[1,1,1,1,1,1],[-5,-4,-3,-2,-1,0],[25/2,16/2,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]]=[[0],[0],[1],[0],[0]]`

on further solving we get the values of the unknown i.e.

Values = a=-0.833 b=5.0833 c=13.00 d=17.833 e=-12.833 f=3.75

These coefficients will further be used in code to get the value of derivative.

close all
clc
%Analytial function defintion and value assignment
x=pi/3;
dx=linspace(pi/10,pi/1000,100);
f=@(x)(exp(x)*cos(x));
g=@(x)(-2*exp(x)*sin(x));
d2f=g(x);
func=f(x);
%For Central Difference Scheme
%Coefficients taken from Taylor's Table
P=[1 1 1 1 1
    -2 -1 0 1 2
    4/2 1/2 0 1/2 4/2
    -8/6 -1/6 0 1/6 8/6
    16/24 1/24 0 1/24 16/24];
C=[0;0;1;0;0];
% As P*a=C i.e a=inverse(P)*c
cf1=P*C;
% For Right Skewed difference scheme
Q=[1 1  1 1 1 1 
    0 1 2 3 4 5 
    0 1/2 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];
D=[0;0;1;0;0;0];
% As Q*b=s i.e b=inverse(Q)*s
cf2=Q*D;

%For Left Skewed difference scheme
R=[1 1 1 1 1 1
    -5 -4 -3 -2 -1 0
    25/2 16/2 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];
E=[0;0;1;0;0;0];
% As R*e=c i.e c=inverse(R)*e
cf3=R*E;

%Calculating values of derivative using CDS RSDC LSDC
%Defining variable size before using in loop
central_diff=zeros(1,100);
err_central=zeros(1,100);
right_sdiff=zeros(1,100);
err_right=zeros(1,100);
left_sdiff=zeros(1,100);
err_left=zeros(1,100);

for i=1:length(dx)
    %Calculating value and error by central difference scheme
    central_diff(i)=(cf1(1)*f(x-2*dx(i))+cf1(2)*f(x-dx(i))+cf1(3)*f(x)+cf1(4)*f(x+dx(i))+cf1(5)*f(x+2*dx(i)))/(dx(i)^2);
    err_central(i)=abs(central_diff(i)-d2f);

    %Calculating value and error by Right skewed difference scheme
    right_sdiff(i)=(cf2(1)*f(x)+cf2(2)*f(x+dx(i))+cf2(3)*f(x+2*dx(i))+cf2(4)*f(x+3*dx(i))+cf2(5)*f(x+4*dx(i))+cf2(6)*f(x+5*dx(i)))/(dx(i)^2);
    err_right(i)=abs(right_sdiff(i)-d2f);

    %Calculating value and error by Left skewed difference scheme
    left_sdiff(i)=(cf3(1)*f(x-5*dx(i))+cf3(2)*f(x-4*dx(i))+cf3(3)*f(x-3*dx(i))+cf3(4)*f(x-2*dx(i))+cf3(5)*f(x-dx(i))+cf3(6)*f(x))/(dx(i)^2);
    err_left(i)=abs(left_sdiff(i)-d2f);
end

%Plotting values of derivative
subplot(2,1,1)
hold on
fplot(d2f,'--')
xlim([-0.05 0.35]);
plot(dx,central_diff,'linewidth',2,'Color','r');
plot(dx,right_sdiff,'Linewidth',2,'Color','b');
plot(dx,left_sdiff,'LineWidth',2,'Color','g');
legend('analytical value at x=pi/3','central diff.','right skewed diff.','left skewed diff.','location','best')
xlabel('Grid siza(dx)');
ylabel('Value');
title('Comparison bw value of derivative calculated by diff, order of approximation')
hold off

%Plotting errors in approximation of derivative
subplot(2,1,2)
loglog(dx,err_central,'linewidth',2,'Color','r');
hold on
xlim([0.002 0.4]);
loglog(dx,err_right,'linewidth',2,'color','b');
loglog(dx,err_left,'linewidth',2,'color','g');
legend('central diff.','right skewed diff.','left skewed diff.','location','best')
xlabel('Grid siza(dx)');
ylabel('absolute erro');
title('Comparioson b/w 4th order approximation method for error in calculating d2f/dx2');

Steps Involved:

1. We assign the value of x and range of dx

2. We define an anonymous function of analytical function and its derivative.

3. Use the matrices from above discussed section and solved it using inversion method i.e.

AX=B ------ X=A-1B

and stored the coefficients in cf1,cf2 and cf3.

4. Assigned the rows and column to the appropriate matrices.

5. Loop the order of approximation expression in a for loop under the range of dx Here the absolute error is calculated from the derivative value calculated analytically.

6. The calculated valued are plotted to do a better comparison.

Result:

The Skewed method is useful for finding the numerical derivatives more accurately. We can see that the central and right error are similar and also there is a hype in the value in the error dx graph. So Skewed schemes are better tha CD method.