Higher Order Approximations using SKEWED AND SYMMETRIC STENCILS
AIM:- Write a MATLAB code to compute 4th order approximations of the Second order derivative of CENTRAL DIFFERENCE, SKEWED RIGHT SIDED and SKEWED LEFT SIDED DIFFERENCE.
THEORY:
Analytical Function - `e^x*cos(x)`
The Second order derivative of the analytical function is `-2*e^(x)*sinx`
CENTRAL DIFFERENCE SCHEME (CDS):
In central differencing, the values are to be taken from both left and right side of the point in the stencil.
So, the Taylor Table for Central difference approximation is
|
`f(x)` |
`trianglex*f'(x)` |
`trianglex^2*f''(x)` |
`trianglex^3*f'''(x)` |
`trianglex^4*f^4(x)` |
|
|
`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 |
4e/2 |
8e/6 |
16e/24 |
|
0 |
0 |
1 |
0 |
0 |
The Taylor series for the central differencing scheme is as follows:
`(d^2u)/(dx^2) = af(i-2) + bf(i-1) + cf(i)+ df(f+1) + ef(i+2)`
From the above table we can solve the equations
`a+b+c+d+e =0`
`-2a-b+d+2e = 0`
`(4a)/2 + b/2+d/2+(4e)/2 = 1`
`(-8a)/6 -b/6+d/6+(8e)/6 =0`
`(16a)/24+b/24+d/24+(16e)/24 = 0`
MATLAB CODE:
%% CENTRAL DIFFERENCE APPROXIMATION
function central_approx = central_diff(x,dx)
A1 = [ 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];
B1 = [0; 0; 1; 0; 0];
C = A1\B1;
% Numerical_Derivative
a = exp(x-2*dx)*cos(x-2*dx);
b = exp(x-dx)*cos(x-dx);
c = exp(x)*cos(x);
d = exp(x+dx)*cos(x+dx);
e = exp(x+2*dx)*cos(x+2*dx);
central_approx = ( (C(1)*a) + (C(2)*b) + (C(3)*c) + (C(4)*d) + (C(5)*e) ) / (dx^2);
end
SKEWED RIGHT DIFFERENCE SCHEME(SRDS):
In skewed right-sided dfference, the values are to be taken only from the right hand side of the point of interest in the stencil.
So, the taylor table for Skewed Right Side is
|
`f(x)` |
`trianglex*f'(x)` |
`trianglex^2*f''(x)` |
`trianglex^3*f'''(x)` |
`trianglex^4*f^4(x)` |
|
|
`af_i` |
a |
0 |
0 |
0 |
0 |
|
`bf_(i+1)` |
b |
b |
b/2 |
b/6 |
b/24 |
|
`cf_(i+2)` |
c |
2c |
4c/2 |
8c/6 |
16c/24 |
|
`df_(i+3)` |
d |
3d |
9d/2 |
27d/6 |
81d/24 |
|
`ef_(i+4)` |
e |
4e |
16e/2 |
64e/6 |
256e/24 |
|
0 |
0 |
1 |
0 |
0 |
The Taylor series for the Skewed right hand side differencing scheme is as follows:
`(d^2u)/(dx^2) = af(i) + bf(i+1) + cf(i+2)+ df(f+3) + ef(i+4)`
From the above table (Skewed right side), we can solve the equations
`a+b+c+d+e =0`
`b+2c+3d+4e= 0`
`b/2+(4c)/2+(9d)/2+(16e)/2 = 1`
`b/6+(8c)/6+(27d)/6+(64e)/6 =0`
`b/24+(16c)/24+(81d)/24+(256e)/24 = 0`
MATLAB CODE:
%% SKEWED RIGHT SIDED DIFFERENCE
function skewed_right_approx = skewed_right(x,dx)
A2 = [1 1 1 1 1;
0 1 2 3 4;
0 1/2 4/2 9/2 16/2;
0 1/6 8/6 27/6 64/6;
0 1/24 16/24 81/24 256/24];
B2 = [0;0;1;0;0];
R = A2\B2;
% Numerical_Derivative
c = exp(x)*cos(x);
d = exp(x+dx)*cos(x+dx);
e = exp(x+2*dx)*cos(x+2*dx);
f = exp(x+3*dx)*cos(x+3*dx);
g = exp(x+4*dx)*cos(x+4*dx);
skewed_right_approx = ( (R(1)*c) + (R(2)*d) + (R(3)*e) + (R(4)*f) + (R(5)*g) ) / (dx^2);
endSKEWED LEFT DIFFERENCE SCHEME(SLDS):
In skewed left-sided dfference, the values are to be taken only from the left hand side of the point of interest in the stencil.
So, the taylor table for Skewed Left Side is
|
`f(x)` |
`trianglex*f'(x)` |
`trianglex^2*f''(x)` |
`trianglex^3*f'''(x)` |
`trianglex^4*f^4(x)` |
|
|
`af_(i-4)` |
a |
-4a |
16a/2 |
-64a/6 |
-256a/24 |
|
`bf_(i-3)` |
b |
-3b |
9b/2 |
-27b/6 |
81b/24 |
|
`cf_(i-2)` |
c |
-2c |
4c/2 |
-8c/6 |
16c/24 |
|
`df_(i-1)` |
d |
-d |
d/2 |
-d/6 |
d/24 |
|
`ef_(i)` |
e |
0 |
0 |
0 |
0 |
|
0 |
0 |
1 |
0 |
0 |
The Taylor series for the Skewed right hand side differencing scheme is as follows:
`(d^2u)/(dx^2) = af(i-4) + bf(i-3) + cf(i-2)+ df(f-1) + ef(i)`
From the above table (Skewed right side), we can solve the equations
`a+b+c+d+e =0`
`-4a-3b-2c-d=0`
`(16a)/2 +(9b)/2+(4c)/2+(d)/2 = 1`
`(-64a)/6-(27b)/6-(8c)/6-d/6=0`
`(256a)/24+(81b)/24+(16c)/24+(d)/24 = 0`
MATLAB CODE:
%% SKEWED LEFT SIDED DIFFERENCE
function skewed_left_approx = skewed_left(x,dx)
A3 = [ 1 1 1 1 1;
-4 -3 -2 -1 0;
16/2 9/2 4/2 1/2 0;
-64/6 -27/6 -8/6 -1/6 0;
256/24 81/24 16/24 1/24 0];
B3 = [0;0;1;0;0];
L = A3\B3;
% Numerical_Derivative
j = exp(x-4*dx)*cos(x-4*dx);
h = exp(x-3*dx)*cos(x-3*dx);
a = exp(x-2*dx)*cos(x-2*dx);
b = exp(x-dx)*cos(x-dx);
c = exp(x)*cos(x);
skewed_left_approx = ( (L(1)*j) + (L(2)*h) + (L(3)*a) + (L(4)*b) + (L(5)*c) ) / (dx^2);
endMAIN CODE:
%% 4th order approximation for 2nd order derivative
% Inputs
x = pi/3;
n = 50;
dx = linspace(pi/4,pi/4000,n);
% Analytical Solution
% Function = exp(x)*cos(x)
Analytical_Solution = -2*exp(x)*sin(x);
central_approx = zeros(1,n);
skewed_right_approx = zeros(1,n);
skewed_left_approx = zeros(1,n);
forward_diff_approx = zeros(1,n);
backward_diff_approx = zeros(1,n);
for i = 1:length(dx)
central_approx(i) = central_diff(x,dx(i));
skewed_right_approx(i) = skewed_right(x,dx(i));
skewed_left_approx(i) = skewed_left(x,dx(i));
forward_diff_approx(i) = forward(x,dx(i));
backward_diff_approx(i) = backward(x,dx(i));
end
% Calculating ERRORS
Central_error = abs(central_approx - Analytical_Solution);
Skewed_right_error = abs(skewed_right_approx - Analytical_Solution);
Skewed_left_error = abs(skewed_left_approx - Analytical_Solution);
Forward_error = abs(forward_diff_approx - Analytical_Solution);
Backward_error = abs(backward_diff_approx - Analytical_Solution);
%Plotting
figure(1),clf
loglog(dx,Central_error,'linewidth',2,'color','r')
hold on
loglog(dx,Skewed_right_error,'linewidth',2,'color','b')
loglog(dx,Skewed_left_error,'linewidth',2,'color','g')
legend('CDS','SRDS','SLDS')
title('Variation of dx w.r.t. CDS,SRDS,SLDS')
xlabel('dx')
ylabel('Error')
%------------------------------------------------------%
figure(2),clf
loglog(dx,Central_error,'linewidth',2,'color','r')
hold on
loglog(dx,Skewed_right_error,'linewidth',2,'color','b')
loglog(dx,Skewed_left_error,'linewidth',2,'color','g')
loglog(dx,Forward_error,'linewidth',2,'color','k')
loglog(dx,Backward_error,'linewidth',2,'color','m')
legend('CDS','SRDS','SLDS','FDS','BDS')
title('Variation of dx w.r.t. CDS,SRDS,SLDS,FDS,BDS')
xlabel('dx')
ylabel('Error')
COMPARISON BETWEEN CDS,SRDS, SLDS
- The graph between dx v/s error shows that the error we get in central difference scheme is very less when compared to other two schemes.
- The reason we get less error is that we consider information from both the sides with reference to the current point, but in other two schemes we take information from only one side of the numerical stencil or only from one direction of the numerical stencil.
------------------------------------------------------------------------------------------------------------
FORWARD DIFFERENCE (FD) :
%% FORWARD DIFFERENCE
function forward_diff_approx = forward(x,dx)
c = exp(x)*cos(x);
d = exp(x+dx)*cos(x+dx);
forward_diff_approx = (d-c)/(dx);
endBACKWARD DIFFERENCE (BD):
%% BACKWARD DIFFERENCE
function backward_diff_approx = backward(x,dx)
c = exp(x)*cos(x);
b = exp(x-dx)*cos(x-dx);
backward_diff_approx = (c-b)/(dx);
end
------------------------------------------------------------------------------------------------------------------
PLOTS AND RESULTS:
COMPARISON BETWEEN CDS, SRDS, SLDS, FD, BD:
- In the graph between error v/s dx, we see the lowest error we get is again from CDS even when compared with FDS and BDS schemes.
- The scale of error from CDS to FD/BD is from 1 to `10^-10` which is very high.
- It concludes that the higher the order of approximation, the lower the error we get (most accurate value we get and closer to the analytical solution).
-------------------------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------------------------
Why a skewed scheme is useful and what can a skewed scheme can do that a CDS cannot do?
Basically, a central difference scheme(CDS) takes a information for current point from the all the neighbouring sides of it ( In X-Direction, it takes information from left and right sides and In Y-Direction, it takes information from top and bottom sides of the current point.)
Where as, a skewed difference takes information for current point from only one side ( for skewed-left it takes information from left side, and for skewed-right it takes information from right side).
So,this helps at boundaries where information only flows in one direction.