Week 3 - Taylor table method and Matlab code

4th order approximations of the second-order derivative:

The given function is exp(x)*cos(x). The same is used to prepare a function so that it can be easily accessed.

function ans = f(x)
  %given function
  ans = exp(x)*cos(x);
endfunction

Central differencing:

This uses 2 points to the left and 2 points to the right side of the point (x) to gather data.

The nodes are calculated using Nodes = p+q-1

where,

p = order of derivative

q = order of approximation

Therefore, Nodes = 2+4-1 = 5

The hypothesis equation:

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

Using Taylor Table Method, a function is prepared to solve a,b,c,d,e values (using matrix multiplication) and substitute the same in the above equation.

function central_difference_value = central_difference_scheme(x,dx)
% Creating matrix from the taylor table
N = [ 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];
%RHS of equation
C = [0;0;1;0;0];
% solving for the constants
A = inv(N)*C;
% main hypothesis equation 
central_difference_value = (A(1,1)*f(x-2*dx)+ A(2,1)*f(x-dx)+ A(3,1)*f(x)+A(4,1)*f(x+dx)+A(5,1)*f(x+2*dx))/(dx^2);
end

Skewed right-sided scheme:

This uses the right side points to the point(x) to gather data.

The nodes are calculated using Nodes = p+q

where,

p = order of derivative

q = order of approximation

Therefore, Nodes = 2+4 = 6

The hypothesis equation is

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

Using Taylor Table Method, a function is prepared to solve a,b,c,d,e,g values (using matrix multiplication) and substitute the same in the above equation.

function right_skewed_value = right_sided_skewed_scheme(x,dx)
% Creating matrix from the taylor table
N = [1 1     1      1       1        1
     0 1     2      3       4        5
     0 1/2   4/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];
%RHS of equation
C = [0;0;1;0;0;0];
% solving for the constants
A = inv(N)*C;
% main hypothesis equation 
right_skewed_value = (A(1,1)*f(x)+ A(2,1)*f(x+dx)+ A(3,1)*f(x+2*dx)+A(4,1)*f(x+3*dx)+A(5,1)*f(x+4*dx)+A(6,1)*f(x+5*dx))/(dx^2);
end

Skewed left-sided scheme:

This uses the left side points to the point(x) to gather data.

The nodes are calculated using Nodes = p+q

where,

p = order of derivative

q = order of approximation

Therefore, Nodes = 2+4 = 6

The hypothesis equation is

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

Using Taylor Table Method, a function is prepared to solve a,b,c,d,e,g values (using matrix multiplication) and substitute the same in the above equation.

function left_skewed_value = left_sided_skewed_scheme(x,dx)
% Creating matrix from the taylor table
N = [ 1         1         1        1       1     1
     -5        -4        -3       -2      -1     0
      25/2      16/2      9/2      4/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 ];
%RHS of equation
C = [0;0;1;0;0;0];
% solving for the constants
A = inv(N)*C;
% main hypothesis equation 
left_skewed_value = (A(1,1)*f(x-5*dx)+ A(2,1)*f(x-4*dx)+ A(3,1)*f(x-3*dx)+A(4,1)*f(x-2*dx)+A(5,1)*f(x-1*dx)+A(6,1)*f(x))/(dx^2);
end

Main Program

Once the values are computed, it should be compared with the analytical value and the same has been plotted.

clear all 
close all
clc
% defining the x and dx
x=pi/3;
dx=linspace(pi/4.5,pi/300,100);
% Value using analytical method
analytical_value = -(2*exp(x)*sin(x))
% calling the user-defined functions and calculating the value
for i = 1 : length(dx)
 left_skewed(i) = left_sided_skewed_scheme(x,dx(i));
 central_difference(i) = central_difference_scheme(x,dx(i));
 right_skewed(i) = right_sided_skewed_scheme(x,dx(i));
end
%Finding the error
for j = 1 : length(dx)
 error_left(j) =  abs(analytical_value - left_skewed(j));
 error_central(j) = abs(analytical_value - central_difference(j));
 error_right(j) =  abs(analytical_value - right_skewed(j));
endfor
%Ploting the graph
figure(1)
loglog(dx,error_left)
hold on 
grid on
loglog(dx,error_right)
loglog(dx,error_central)
hold off 
title('Comparison between central,left and right skewed')
xlabel ('dx')
ylabel('error')
legend('Left side skewed','right side skewed','Central difference',"location","northwest")
axis([pi/300 pi/4.5 0 1])

Result:

The graph:

The above graph 1 shows the error for dx range between pi/300 to pi/4.5.

As shown in graph 1, the error significantly reduces as the dx value reduces.

.

The above graph 2 shows the error for dx range between pi/300000 to pi/4.5. Thus we can say a large reduction in the range of dx will increase the error after one point.

Also, the error due to central differencing is less than right and left-skewed scheme.

Even though the error value of central differencing is low, the skewed scheme is useful when the x point which we are calculating doesn't have left or right points to gather data.

Note:

Make sure the function name and the file name matches.