Week 3 - Taylor table method and Matlab code

Problem Statement:

Derive the following 4th order approximations of the second-order derivative.

1. Central difference

2. Skewed right-sided difference

3. Skewed left-sided difference

After deriving these schemes, prove that your skewed schemes are fourth-order accurate i.e the thumb rule given in HINT (ii).

Once you have done this, write a program in Matlab to evaluate the second-order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations that you have derived.

Provide a plot that compares the absolute error between the above-mentioned schemes

In your own words, describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot?

Theory:

1.Central Difference Scheme: In Central Difference Scheme if we need to find the value of a derivative at a particular point, we consider that particular point and the neighbouring points besides that particular

the point in both sides.If the order of the derivative is 'n' and the order of approximation is 'm' then the number of points to be used is m+n-1.Here in this problem order of the derivative is '2' and the order of

approximation is '4' then the no. of points are 4+2-1 = 5 i.e `f_(i-2),f_(i-1),f_(i),f_(i+1),f_(i+2)`.Now we write a Hypothesis equation with five coefficients a,b,c,d,e, and try to determine the numerical scheme

for the second-order derivative which is of fourth-order approximation.

Hypothesis Equation.

`(d^2u)/(dx^2) = a*f(i-2) +b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2)`

Now we apply Taylor's series for each term in R.H.S side.

`a*f(i-2 ) = a*f(i) - a*f'(i)*((2*dx))/(1!)+a*f''(i)*((2*dx)^2)/(2!)-a*f'''(i)*((2*dx)^3)/(3!)+a*f''''(i)*((2*dx)^4)/(4!)-a*f'''''(i)*((2*dx)^5)/(5!)......`

`b*f(i-1 ) = b*f(i) - b*f'(i)*((dx))/(1!)+b*f''(i)*((dx)^2)/(2!)-b*f'''(i)*((dx)^3)/(3!)+b*f''''(i)*((dx)^4)/(4!)-b*f'''''(i)*((dx)^5)/(5!)......`

`c*f(i) = c*f(i)`

`d*f(i+1 ) = d*f(i) + d*f'(i)*((dx))/(1!)+d*f''(i)*((dx)^2)/(2!)+d*f'''(i)*((dx)^3)/(3!)+d*f''''(i)*((dx)^4)/(4!)+d*f'''''(i)*((dx)^5)/(5!)....`

`e*f(i+2 ) = e*f(i) + e*f'(i)*((2*dx))/(1!)+e*f''(i)*((2*dx)^2)/(2!)+e*f'''(i)*((2*dx)^3)/(3!)+e*f''''(i)*((2*dx)^4)/(4!)+e*f'''''(i)*((2*dx)^5)/(5!)......`

`TAYLOR TABLE`	`f(i)`	`f'(i)`	`f''(i)`	`f'''(i)`	`f''''(i)`	`f'''''(i)`	`f''''''(i)`
`a*f(i-2)`	`a`	`-a*2*dx`	`(a*(2*dx)^2)/(2!)`	`-(a*(2*dx)^3)/(3!)`	`(a*(2*dx)^4)/(4!)`	`-a*f'''''(i)*((2*dx)^5)/(5!)`	`a*f''''''(i)*((2*dx)^6)/(6!)`
`b*f(i-1)`	`b`	`-b*dx`	`(b*(dx)^2)/(2!)`	`-(b*(dx)^3)/(3!)`	`(b*(dx)^4)/(4!)`	`-b*f'''''(i)*((dx)^5)/(5!)`	`b*f''''''(i)*((dx)^6)/(6!)`
`c*f(i)`	`c`	`0`	`0`	`0`	`0`	`0`	`0`
`d*f(i+1)`	`d`	`d*dx`	`(d*(dx)^2)/(2!)`	`(d*(dx)^3)/(3!)`	`(d*(dx)^4)/(4!)`	`d*f'''''(i)*((dx)^5)/(5!)`	`d*f''''''(i)*((dx)^6)/(6!)`
`e*f(i+2)`	`e`	`e*2*dx`	`(e*(2*dx)^2)/(2!)`	`(e*(2*dx)^3)/(3!)`	`(e*(2*dx)^4)/(4!)`	`e*f'''''(i)*((2*dx)^5)/(5!)`	`e*f''''''(i)*((2*dx)^6)/(6!)`
`sum(ALL TERMS)`	`0`	`0`	`1`	`0`	`0`	`?`	`?`

1. `a+b+c+d+e = 0`
2. `-a*2 -b+d+e*2 = 0`
3. `((4/2a)+(1/2b)+(1/2d)+(4/2e) )*dx^2= 1`
4. `-8a-b+d+8e = 0`
5. `16a+b+d+16e = 0`
6. `(-32a-b+d+32e)/120*dx^5`
7. `(64a+b+d+64e)/720*dx^6`

The last equations are just coefficients of the error terms and not necessarily equal to zero.

On solving the above four equations we get

`a = -1/(12*dx^2)`

`b = 16/(12*dx^2)`

`c = -30/(12*dx^2)`

`d = 16/(12*dx^2)`

`e = -1/(12*dx^2)`

When we substitute the above a,bc,d,e values in the 6th equation we get coefficient value equal to zero,and when substituted in equation 7th the coefficient we get -1/90*dx^4.since the coefficient is equal to non zero so the order of accurate is equal to `(dx)^4`.

After substituting all values i.en a,b,c,d,e in the Hypothesis equation we get

`(d^2u)/dx^2 = (-1*f(i-2)+16*f(i-1)-30*f(i)+16*f(i+1)-1*f(i+2))/12*dx^2`

2. Skewed right-sided difference: In Skewed right-sided difference if we need to find the value of a derivative at a particular point, we consider that particular point and the neighbouring points which are right

to that particular point.If the order of the derivative is 'n' and the order of approximation is 'm' then the number of points to be used is m+n.Here in this problem order of the derivative is '2' and the

order of approximation is '4' then the no. of points are 4+2 = 6 i.e `f(i),f(i+1),f(i+2),f(i+3),f(i+4),f(i+5)`.Now we write a Hypothesis equation with five coefficients a,b,c,d,e, and try to determine the numerical

scheme for the second-order derivative which is of fourth-order approximation.

Hypothesis Equation:

`(d^2u)/(dx^2) = a*f(i) +b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+g*f(i+5)`

Now we apply Taylor's series for each term in R.H.S side.

1. `a*f(i) = a*f(i)`
2. `b*f(i+1 ) = b*f(i) + b*f'(i)*((dx))/(1!)+b*f''(i)*((dx)^2)/(2!)+b*f'''(i)*((dx)^3)/(3!)+b*f''''(i)*((dx)^4)/(4!)+b*f'''''(i)*((dx)^5)/(5!)....`
3. `c*f(i+2 ) = c*f(i) + c*f'(i)*((2*dx))/(1!)+c*f''(i)*((2*dx)^2)/(2!)+c*f'''(i)*((2*dx)^3)/(3!)+c*f''''(i)*((2*dx)^4)/(4!)+c*f'''''(i)*((2*dx)^5)/(5!)......`
4. `d*f(i+3 ) = d*f(i) + d*f'(i)*((3*dx))/(1!)+d*f''(i)*((3*dx)^2)/(2!)+d*f'''(i)*((3*dx)^3)/(3!)+d*f''''(i)*((3*dx)^4)/(4!)+d*f'''''(i)*((3*dx)^5)/(5!)......`
5. `e*f(i+4 ) = e*f(i) + e*f'(i)*((4*dx))/(1!)+e*f''(i)*((4*dx)^2)/(2!)+e*f'''(i)*((4*dx)^3)/(3!)+e*f''''(i)*((4*dx)^4)/(4!)+e*f'''''(i)*((4*dx)^5)/(5!)......`
6. `g*f(i+5 ) = g*f(i) + g*f'(i)*((5*dx))/(1!)+g*f''(i)*((5*dx)^2)/(2!)+g*f'''(i)*((5*dx)^3)/(3!)+g*f''''(i)*((5*dx)^4)/(4!)+g*f'''''(i)*((5*dx)^5)/(5!)......`

	`f(i)`	`f'(i)`	`f''(i)`	`f'''(i)`	`f''''(i)`	`f'''''(i)`
`a*f(i)`	`a`	`0`	`0`	`0`	`0`	`0`
`b*f(i+1)`	`b`	`b`	`1/2*b`	`b`	`b`	`b`
`c*f(i+2)`	`c`	`c*2`	`4/2*c`	`8*c`	`16*c`	`32*c`
`d*f(i+3)`	`d`	`d*3`	`9/2*d`	`27*d`	`81*d`	`243*d`
`e*f(i+4)`	`e`	`e*4`	`16/2*e`	`64*e`	`256*e`	`1024*e`
`g*f(i+5)`	`g`	`5*g`	`25/2*g`	`125*g`	`625*g`	`3125*g`
`sum`	`0`	`0`	`1`	`0`	`0`	`0`

`a+b+c+d+e+g = 0`

`b+2*c+3*d+4*e+5*g = 0`

`((1/2)*b+(4/2)*c+(9/2)*d+(16/2)*e+(25/2)*g)*dx^2 = 1`

`b+(8)*c+(27)*d+(64)*e+125*g = 0`

`b+16*c+81*d+256*e+625*g = 0`

`b+32*c+243*d+1024*e+3125*g = 0`

On solving above equations:

we get `a = (45/(12*dx^2))`

`b = (-154/(12*dx^2))`

`c = (214/(12*dx^2))`

`d = (-156/(12*dx^2))`

`e = (61/(12*dx^2))`

`g = (-10/(12*dx^2))`

After substituting all values i.en a,b,c,d,e in the Hypothesis equation we get

`(d^2u)/dx^2 = (45*f(i)-154*f(i+1)+214*f(i+2)-156*f(i+3)+61*f(i+4)-10*f(i+5))/(12*dx^2)`

Now we to find the order of accuracy and to do that we have to find the coefficient of `f''''''`term and the term is `(b+64*c+729*d+4096*e+15625*g)/720*dx^6`After substituting all values in the `f''''''`term

we get the coefficient value as `-548*dx^4`.Since the coefficient is non-zero and the power od dx term is 4, so the scheme is fourth-order accurate.

3. Skewed left-sided difference:

In Skewed left-sided difference if we need to find the value of a derivative at a particular point, we consider that particular point and the neighbouring points which are left to

that particular point.If the order of the derivative is 'n' and the order of approximation is 'm' then the number of points to be used is m+n.Here in this problem order of the derivative is '2' and the

order of approximation is '4' then the no. of points are 4+2 = 6 i.e `f(i),f(i-1),f(i-2),f(i-3),f(i-4),f(i-5)`.Now we write a Hypothesis equation with five coefficients a,b,c,d,e, and try to determine the numerical

scheme for the second-order derivative which is of fourth-order approximation.

Hypothesis Equation:

`(d^2u)/(dx^2) = a*f(i) +b*f(i-1)+c*f(i-2)+d*f(i-3)+e*f(i-4)+g*f(i-5)`

Now we apply Taylor's series for each term in R.H.S side.

`a*f(i) = a*f(i)`

`b*f(i-1 ) = b*f(i) - b*f'(i)*((dx))/(1!)+b*f''(i)*((dx)^2)/(2!)-b*f'''(i)*((dx)^3)/(3!)+b*f''''(i)*((dx)^4)/(4!)-b*f'''''(i)*((dx)^5)/(5!)......`

`c*f(i-2 ) = c*f(i) - c*f'(i)*((2*dx))/(1!)+c*f''(i)*((2*dx)^2)/(2!)-c*f'''(i)*((2*dx)^3)/(3!)+c*f''''(i)*((2*dx)^4)/(4!)-c*f'''''(i)*((2*dx)^5)/(5!)......`

`d*f(i-3 ) = d*f(i) - d*f'(i)*((3*dx))/(1!)+d*f''(i)*((3*dx)^2)/(2!)-d*f'''(i)*((3*dx)^3)/(3!)+d*f''''(i)*((3*dx)^4)/(4!)-d*f'''''(i)*((3*dx)^5)/(5!)......`

`e*f(i-4 ) = e*f(i) - e*f'(i)*((4*dx))/(1!)+e*f''(i)*((4*dx)^2)/(2!)-e*f'''(i)*((4*dx)^3)/(3!)+e*f''''(i)*((4*dx)^4)/(4!)-e*f'''''(i)*((4*dx)^5)/(5!)......`

`g*f(i-5 ) = g*f(i) - g*f'(i)*((5*dx))/(1!)+g*f''(i)*((5*dx)^2)/(2!)-g*f'''(i)*((5*dx)^3)/(3!)+g*f''''(i)*((5*dx)^4)/(4!)-g*f'''''(i)*((5*dx)^5)/(5!)......`

	`f(i)`	`f'(i)`	`f''(i)`	`f'''(i)`	`f''''(i)`	`f'''''(i)`
`a*f(i)`	`a`	`0`	`0`	`0`	`0`	`0`
`b*f(i-1)`	`b`	`-b`	`(1/2)*b`	`-b`	`b`
`c*f(i-2)`	`c`	`-2*c`	`(4/2)*c`	`-8*c`	`16*c`
`d*f(i-3)`	`d`	`-3*d`	`(9/2)*d`	`-27*d`	`81*d`
`e*f(i-4)`	`e`	`-4*e`	`(16/2)*e`	`-64*e`	`256*e`
`g*f(i-5)`	`g`	`-5*g`	`(25/2)*g`	`125*g`	`625*g`
	`0`	`0`	`1`	`0`	`0`

`a+b+c+d+e+g = 0`

`-(b+2*c+3*d+4*e+5*g )= 0`

`((1/2)*b+(4/2)*c+(9/2)*d+(16/2)*e+(25/2)*g)*dx^2 = 1`

`-(b+(8)*c+(27)*d+(64)*e+125*g )= 0`

`b+16*c+81*d+256*e+625*g = 0`

`-(b+32*c+243*d+1024*e+3125*g) = 0`

On solving above equation we get

`a = (45/(12*dx^2))`

`b = (-154/(12*dx^2))`

`c = (214/(12*dx^2))`

`d = (-156/(12*dx^2))`

`e = (61/(12*dx^2))`

`g = (-10/(12*dx^2))`

After substituting all values i.en a,b,c,d,e in the Hypothesis equation we get

`(d^2u)/dx^2 = (45*f(i)-154*f(i-1)+214*f(i-2)-156*f(i-3)+61*f(i-4)-10*f(i-5))/(12*dx^2)`

Now we to find the order of accuracy and to do that we have to find the coefficient of `f''''''`term and the term is `(b+64*c+729*d+4096*e+15625*g)/720*dx^6`After substituting all values in the `f''''''`term

we get the coefficient value as `-548*dx^4`.Since the coefficient is non-zero and the power of dx term is 4, so the scheme is fourth-order accurate.

1. Central difference Function Matlab Program:

function fourth_order_cds_error = fourth_order_cds(x,dx)

%%Ananlytical Function = exp(x)*cos(x)

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

 fourth_order_cds = (-1*cos(x-(2*dx))*exp(x-(2*dx)))+16*cos(x-(1*dx))*exp(x-(1*dx))-30*cos(x)*exp(x)+16*cos(x+(1*dx))*exp(x+(1*dx))-1*cos(x+(2*dx))*exp(x+(2*dx))/12*dx^2

 fourth_order_cds_error = abs( fourth_order_cds-Ananlytical_derivative);
 
end

2. Skewed right-sided difference Function Matlab Program:

function fourth_order_srs_error = fourth_order_srs(x,dx)

%%Ananlytical Function = exp(x)*cos(x)

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

 fourth_order_srs = ((45*cos(x*dx)*exp(x*dx))-(154*cos(x+(1*dx))*exp(x+(1*dx)))+(214*cos(x+(2*dx))*exp(x+(2*dx)))-(156*cos(x+(3*dx))*exp(x+(3*dx)))+(61*cos(x+(4*dx))*exp(x+(4*dx)))-(10*cos(x+(5*dx))*exp(x+(5*dx))))/12*dx^2
 
 fourth_order_srs_error = abs(fourth_order_srs-Ananlytical_derivative);
 
end

3. Skewed left-sided difference Function Matlab Program:

function fourth_order_sls_error = fourth_order_sls(x,dx)

%%Ananlytical Function = exp(x)*cos(x)

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

 fourth_order_sls  = ((45*cos(x*dx)*exp(x*dx))-(154*cos(x-(1*dx))*exp(x-(1*dx)))+(214*cos(x-(2*dx))*exp(x-(2*dx)))-(156*cos(x-(3*dx))*exp(x-(3*dx)))+(61*cos(x-(4*dx))*exp(x-(4*dx)))-(10*cos(x-(5*dx))*exp(x-(5*dx))))/(12*dx^2);

 fourth_order_sls_error = abs(fourth_order_sls-Ananlytical_derivative);
 
end

Main MATLAB Program:

%%write a program in Matlab to evaluate the second-order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations that you have derived.

clear all %clear command window
close all %close workspace
clc

x = pi/3;%Intial point
dx= linspace(pi/10,pi/3000,200)%Defining the domain for dx.

%%Ananlytical Function = exp(x)*cos(x)
Ananlytical_derivative = -2*exp(x)*sin(x)

for i = 1:length(dx)
    
     fourth_order_cds_error (i) = fourth_order_cds(x,dx(i));
     fourth_order_srs_error(i) = fourth_order_srs(x,dx(i));
     fourth_order_sls_error(i) = fourth_order_sls(x,dx(i));
    
end

figure(1)
loglog(dx,fourth_order_cds_error,'g')%logarithmic scales are used for both the X- and Y- axes.
hold on
loglog(dx,fourth_order_srs_error,'b')%logarithmic scales are used for both the X- and Y- axes.
hold on
loglog(dx,fourth_order_sls_error,'r')%logarithmic scales are used for both the X- and Y- axes.
title('comaparing dx and Absolute Error')%adds the specified title
legend('central scheme','skewed Right scheme','skewed left scheme')%creates a legend with descriptive labels
xlabel('dx')%adds text beside the X-axis
ylabel('Absolute Error')%adds text beside the y-axis

Output:

Step-by-step procedure:

Based upon the type of scheme define the neighbouring points and define the Hypothesis equation using variables corresponding to the neighbouring points. Expand each term in the hypothesis equation with

Taylor's series and add corresponding derivatives of all terms. The number of equations must be equal to the number of unknown coefficients and sum it up to zero. The equations are solved to get values of all

coefficients.

1.A function is created named as fourth_order_cds which calculates the fourth-order approximation using central differencing scheme

2.A function is created named as fourth_order_sls which calculates the fourth-order approximation using skewed left the differencing scheme.

3.A function is created named as fourth_order_srs which calculator the fourth-order approximation using skewed right the differencing scheme.

4.All the three functions are called into the main program. And define the current point 'x' and space 'dx' which specifies the distance between neighbouring points.

5.The error between the analytical solution and approximate solution by using Matlab command ‘abs’.

6. A graph is drawn between the error and order of approximation by using Matlab command ‘log-log’.

Conclusion:

1.Through central difference scheme we get a low error then the other two schemes and in the central difference scheme, the information is used from both sides.

2.Although CDS gives more accurate approximation but it cannot simulate at boundaries because it requires information from both sides.

3.The skewed schemes can operate in boundary as it can operate inflow of information either in right or left side.

4. As the 'dx' value increases the error increases up to a point and then decreases due to roundoff error.