Taylor table method and Matlab code to calculate 4th order approximations of 2nd derivative

Objective:

• Derive the following 4th order approximations of the second-order derivative.
  1. Central difference
  2. Skewed right-sided difference
  3. Skewed left-sided difference
• Also, we have to prove that the skewed schemes are fourth-order accurate.
  
  
• Write a program in Matlab to evaluate the second-order derivative of the analytical function `f(x) = e^x*cos(x)` and compare it with the 3 numerical approximations that we have derived.
• Using a plot, compare the absolute error with respect to the above-mentioned schemes, and also explain the usefulness of the skewed scheme over the Central Differencing scheme.

Solution:

The derivative of a function "f" at the point "x" is defined as the limit of a difference quotient:

`f'(x) = lim_(hto0) (f(x+h) - f(x))/h`

In other words, the difference quotient `(f(x+h) - f(x))/h` is an approximation of the derivative `f'(x)`, and this approximation gets better as `h` gets smaller.

Taylor’s theorem with remainder gives the Taylor series expansion:

`f(x+h) = f(x) + h*f'(x) + h^2*(f''(x+ Deltah))/(2!) +...`

Rearranging this gives us:

`f'(x)= (f(x+h) - f(x))/h + h*f''(x+h)/(2)`

which tells us that the error is proportional to h to the power 1, so `(f(x+h) - f(x))/h` is said to be a “first-order” approximation.

To derive a higher-order approximation, we need to Discretize the PDE in order to solve it on MATLAB. For this, we are going to use Taylor's Table Method to convert the PDE in a linear equation.

The Taylor's Table is just a tabulation of expansion of every point of approximation taken, this helps us to find the exact values of coefficients that are required to make the PDE a linear equation.

To form Taylor's Table we first need a stencil to work on, and it will depend on the type of approximation we are taking. It can be a central approximation where the value at a given point is calculated using the points on its both side or it can be biased or skewed to any side.

Now, in order to calculate the error in each type of approximation, we first need to form Taylor's Table for each type, then we have to calculate the values of coefficients and then find the approximate value which is then compared to the analytical value.

The basic stencil for Central differencing, Left-skewed & Right-skewed Taylor table is shown below:

• 4th Order Central Differencing approximation of 2nd derivative of a function:

• Based on the stencil, the approximate equation will be:

`(del^2f)/(delx^2)approx`

`"{ "u*exp(x-2*dx)*cos(x-2*dx) + v*exp(x-dx)*cos(x-dx) + g*exp(x)*cos(x) +`

`(y*exp(x+dx)*cos(x+dx) + z*exp(x+2*dx)*cos(x+2*dx)" }")/(dx^2)`

Taylor's table for the central difference stencil to calculate 4th order approximation of 2nd derivative of a given function will be:

The summation of all equations will give an approximate answer to the 2nd derivative of the function.

we have 5 unknows and to compute them we need 5 equations. From Taylor's Table, the equation formed from the `f''(x)` will be equal to 1 as we are calculating approximation for the 2nd order derivative and therefor, other 4 equations except for `f^v(x) " & " f^(vi)(x)` will be equal to zero.

The equations formed from Taylor's Table are:

1. `u+v+g+y+z = 0`
2. `-2*u-v+0+y+2*z=0`
3. `4/2*u+1/2*v+0+1/2*y+4/2*z=1`
4. `-8/6*u-1/6*v+0+1/6*y+8/6*z=0`
5. `16/24*u+1/24*v+0+1/24*y+16/24*z=0`
6. `-32/120*u-1/120*v+0+1/120*y+32/120*z=?`
7. `64/240*u+1/240*v+0+1/240*y+64/240*z=?`

Matrix for central difference Taylor table is given below (ignoring equation VI & VII):

`"central_difference_m" = [[1,1,1,1,1],[-2,-1,0,1,2],[2,1/2,0,1/2,2],[-8/6,-1/6,0,1/6,8/6],[16/24,1/24,0,1/24,16/24]]`

`"b_c" = [[0],[0],[1],[0],[0]]`

On solving this matrix in MATLAB we get the following results for the coefficients:

central_diffrence_coeff=inv(central_diffrence_m)*b_c;

u= -0.0833

v= 1.3333

g= -2.5000

y= 1.3333

z= -0.0833

On analysing the value of the coefficients, we can see that `u=z " & "v=y` and on using these values in equation VI, we can say that `-32/120*u-1/120*v+0+1/120*y+32/120*z=0`.

Now, on rewriting the equation we get:

`(del^2f)/(delx^2)=f''(x)*Deltax^2 + (64/240*u+1/240*v+0+1/240*y+64/240*z)*f^(vi)*Deltax^6`

which after taking `Deltax^2` common, can be rewritten as:

`(del^2f)/(delx^2)approx`

`"{ "u*exp(x-2*dx)*cos(x-2*dx) + v*exp(x-dx)*cos(x-dx) + g*exp(x)*cos(x) +`

`(y*exp(x+dx)*cos(x+dx) + z*exp(x+2*dx)*cos(x+2*dx)" }")/(dx^2)`

`= f''(x)+ O.Deltax^4`

Now that we are getting the Truncation error term of 4th order, it can be said that this is a 4th order approximation of 2nd derivative of a function.

• 4th Order Right-Skewed approximation of 2nd derivative of a function:

• Based on the stencil, the approximate equation will be:

`(del^2f)/(delx^2)approx`

`"{ "(a*exp(x)*cos(x) + b*exp(x+dx)*cos(x+dx) + c*exp(x+2*dx)*cos(x+2*dx) +`

`(d*exp(x+3*dx)*cos(x+3*dx) + e*exp(x+4*dx)*cos(x+4*dx) + h*exp(x+5*dx)*cos(x+5*dx) " }")/(dx^2)`

Taylor's table for the right-skewed stencil to calculate 4th order approximation of 2nd derivative will be:

The summation of all equations will give an approximate answer to the 2nd derivative of the function.

we have 6 unknows and to compute them we need 6 equations. From Taylor's Table, the equation formed from the `f''(x)` will be equal to 1 as we are calculating approximation for the 2nd order derivative and therefor, other 5 equations except for `f^(vi)(x)` will be equal to zero.

The equations formed from Taylor's Table are:

1. `a+b+c+d+e+h=0`
2. `0+b + 2*c + 3*d + 4*e + 5*h = 0`
3. `0+b/2 + 2*c + 9*d/2 + 16*e/2 + 25*h/2 = 1`
4. `0+b/6 + 8*c/6 + 27*d/6 + 64*e/6 + 125*h/6 = 0`
5. `0+b/24 + 16*c/24 + 81*d/24 +256*e/24 + 625*h/24 = 0`
6. `0+b/120 + 32*c/120 + 243*d/120 + 1024*e/120 + 3125*h/120 = 0`
7. `0+b/720+64*c/720+729*d/720+4096*e/720+15625*h/720=?`

Matrix for right-skewed Taylor table is given below (ignoring equation VII):

`"right_skewed_m" = [[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]]`

`b=[[0],[0],[1],[0],[0],[0]]`

On solving this matrix in MATLAB we get the following results for the coefficients:

right_skewed_coeff=inv(right_skewed_m)*b;

a= 3.7500

b= -12.8333

c= 17.8333

d= -13.0000

e= 5.0833

h= -0.8333

Now, on rewriting the equation we get:

`(del^2f)/(delx^2)=f''(x)*Deltax^2 + (b/720+64*c/720+729*d/720+4096*e/720+15625*h/720)*f^(vi)*Deltax^6`

which after taking `Deltax^2` common, can be rewritten as:

`(del^2f)/(delx^2)approx`

`"{ "a*exp(x)*cos(x) + b*exp(x+dx)*cos(x+dx) + c*exp(x+2*dx)*cos(x+2*dx) +`

`(d*exp(x+3*dx)*cos(x+3*dx) + e*exp(x+4*dx)*cos(x+4*dx) + h*exp(x+5*dx)*cos(x+5*dx)" }")/(dx^2)`

`= f''(x)+ O.Deltax^4`

Now that we are getting the Truncation error term of 4th order, it can be said that this is a 4th order approximation of 2nd derivative of the function.

• 4th Order Right-Skewed approximation of 2nd derivative of a function:

• Based on the stencil, the approximate equation will be:

`(del^2f)/(delx^2)approx`

`"{ "m*exp(x-5*dx)*cos(x-5*dx) + n*exp(x-4*dx)*cos(x-4*dx) + o*exp(x-3*dx)*cos(x-3*dx) +`

`(p*exp(x-2*dx)*cos(x-2*dx) + q*exp(x-1*dx)*cos(x-1*dx) + r*exp(x)*cos(x)" }")/(dx^2)`

Taylor's table for the left-skewed stencil to calculate 4th order approximation of 2nd derivative will be:

The summation of all equations will give an approximate answer to the 2nd derivative of the function.

we have 6 unknows and to compute them we need 6 equations. From Taylor's Table, the equation formed from the `f''(x)` will be equal to 1 as we are calculating approximation for the 2nd order derivative and therefor, other 5 equations except for `f^(vi)(x)` will be equal to zero.

The equations formed from Taylor's Table are:

1. `m+n+"o"+p+q+r=0`
2. `-5*m-4*n-3*o-2*q-0=0`
3. `25*m/2+16*n/2+9*o/2+2*p+q/2+0=0`
4. `-125*m/6-64*n/6-27*o/6-8*p/6-q/6-0=0`
5. `625*m/24+256*n/24+81*o/24+16*p/24+q/24+0=0`
6. `-3125*m/120-1024*n/120-243*o/120-32*p/120-q/120-0=0`
7. `15625*m/720+4096*n/720+729*o/720+64*p/720+q/720+0=0`

Matrix for left-skewed Taylor table is given below (ignoring equation VII):

`"left_skewed_m" = [[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]]`

On solving this matrix in MATLAB we get the following results for the coefficients:

left_skewed_coeff=inv(left_skewed_m)*b;

m= -0.8333

n= 5.0833

o= -13.0000

p= 17.8333

q= -12.8333

r= 3.7500

Now, on rewriting the equation we get:

`(del^2f)/(delx^2)=f''(x)*Deltax^2 + (15625*m/720+4096*n/720+729*o/720+64*p/720+q/720)*f^(vi)*Deltax^6`

which after taking `Deltax^2` common, can be rewritten as:

`(del^2f)/(delx^2)approx`

`"{ "m*exp(x-5*dx)*cos(x-5*dx) + n*exp(x-4*dx)*cos(x-4*dx) + o*exp(x-3*dx)*cos(x-3*dx) +`

`(p*exp(x-2*dx)*cos(x-2*dx) + q*exp(x-1*dx)*cos(x-1*dx) + r*exp(x)*cos(x)" }")/(dx^2)`

`= f''(x)+ O.Deltax^4`

Now that we are getting the Truncation error term of 4th order, it can be said that this is a 4th order approximation of 2nd derivative of the function.

• MATLAB programming:

Function: `f(x) = e^x*cos(x)`

at `x=pi/3`

where Dx is calculated by the following command:

Dx=linspace(pi/400,pi/15,50);

This means that we are taking 50 values in-between `pi/400 " & " pi/15`

In the MATLAB we defined the matrix and calculated the values of coefficient and directly used them in our approximation formulas in order to remove any approximation errors that might affect the values of the coefficient which in turn will produce more error in the final answer.

The main program to calculate the analytical derivative, Left-skewed, Right-skewed, Central Differencing of 4th order and then calculate the absolute error and then plot the graph between absolute error and Dx is given below:

clear all
close all
clc
%calculating right_skewed_coefficients
right_skewed_m=[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];
B=[0;0;1;0;0;0];
right_skewed_coeff=inv(right_skewed_m)*B;
a=right_skewed_coeff(1)
b=right_skewed_coeff(2)
c=right_skewed_coeff(3)
d=right_skewed_coeff(4)
e=right_skewed_coeff(5)
h=right_skewed_coeff(6)

% calculating left_skewed_coefficients
left_skewed_m = [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];
left_skewed_coeff = inv(left_skewed_m)*B;
m=left_skewed_coeff(1)
n=left_skewed_coeff(2)
o=left_skewed_coeff(3)
p=left_skewed_coeff(4)
q=left_skewed_coeff(5)
r=left_skewed_coeff(6)

% calculating central_diffrence_coefficients
central_diffrence_m=[1,1,1,1,1;-2,-1,0,1,2;2,1/2,0,1/2,2;-8/6,-1/6,0,1/6,8/6;16/24,1/24,0,1/24,16/24];
b_c=[0;0;1;0;0];
central_diffrence_coeff=inv(central_diffrence_m)*b_c;
u=central_diffrence_coeff(1)
v=central_diffrence_coeff(2)
g=central_diffrence_coeff(3)
y=central_diffrence_coeff(4)
z=central_diffrence_coeff(5)

% defining function, value of x and dx
%f(x) = (exp(x))*cos(x)
x=pi/3;

Dx=linspace(pi/400,pi/15,50);

%evaluating the 4th order approximation of 2nd derivative of function using for loop

for i=1:length(Dx)
dx=(Dx(i));
analytical_second_derivative(i)=((-2)*exp(x))*sin(x);
left_skewed_approx(i) = (m*exp(x-5*dx)*cos(x-5*dx) + n*exp(x-4*dx)*cos(x-4*dx) + o*exp(x-3*dx)*cos(x-3*dx) + p*exp(x-2*dx)*cos(x-2*dx) + q*exp(x-1*dx)*cos(x-1*dx) + r*exp(x)*cos(x))/(dx^2);
left_skewed_error(i) = abs(left_skewed_approx(i)-analytical_second_derivative(i));
right_skewed_approx(i) = (a*exp(x)*cos(x) + b*exp(x+dx)*cos(x+dx) + c*exp(x+(2*dx))*cos(x+(2*dx)) + d*exp(x+(3*dx))*cos(x+(3*dx)) + e*exp(x+(4*dx))*cos(x+(4*dx)) + h*exp(x+(5*dx))*cos(x+(5*dx)))/(dx^2);
right_skewed_error(i) = abs(right_skewed_approx(i)-analytical_second_derivative(i));
central_differences_approx(i) = ((u*exp(x-(2*dx))*cos(x-(2*dx))) + (v*exp(x-dx)*cos(x-dx)) + (g*exp(x)*cos(x)) + (y*exp(x+dx)*cos(x+dx)) + (z*exp(x+(2*dx))*cos(x+(2*dx))))/(dx^2);
central_error(i) = abs(central_differences_approx(i)-analytical_second_derivative(i));
end
analytical_second_derivative

%ploting the calculated error vs Dx graph
figure(1)
plot(Dx,left_skewed_error,'-*',Dx,right_skewed_error,'-p',Dx,central_error,'-s')
title(['plot of 2nd derivative approximate of F(x) = exp(x)*cos(x), at x=',num2str(x)])
xlabel('value of dx')
ylabel('Absolute error compared to analytical value')
legend('left skewed error','right skewed error','central error')

%ploting log graph of calculated error vs Dx to get a better understanding
figure(2)
loglog(Dx,left_skewed_error,'-*',Dx,right_skewed_error,'-p',Dx,central_error,'-s')
title(['plot of 2nd derivative approximate of F(x) = exp(x)*cos(x), at x=',num2str(x)])
xlabel('value of dx')
ylabel('Absolute error compared to analytical value')
legend('left skewed error','right skewed error','central error')

In the above program, we used a for loop to calculate the value of approximations and absolute error compared to the analytical solution at a different value of Dx.

A normal plot is used to visually compare the errors of all 3 types of approximations, but this doesn't tell much about the error in the Central Differencing scheme, therefore a logarithmic pot is used to understand the errors more.

Comparison of absolute error of 4th order approximation of 2nd derivative of the function is shown below:

A logarithmic comparison of absolute error of 4th order approximation of 2nd derivative of the function to get a better understanding of absolute error is shown below:

From the above plots, we can clearly say that the value of absolute error will decrease as the value of Dx decrease in all schemes, but it can be noted that the Central Difference scheme start giving error of the order of `10^-4` from the very beginning, whereas the other 2 schemes get this order of error when the value of Dx gets nearly less than half of that in the first case of CDS.

This clearly shows the advantage of using CDS approximation over the Left-Skewed or Right-Skewed scheme.

Ques.: Why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?

Ans.: Even though we are getting much more error in Left-Skewed & Right-Skewed approximation scheme, it is very useful in the condition of the boundary of the domain where there are not enough points available to apply Central Differencing approximation scheme as it requires grid points on both sides of the reference point. In such cases, only Left-Skewed or Right-Skewed approximation scheme can be applied.