Week 3 - Taylor table method and Matlab code
1. AIM:
To derive the 4th order approximation of the second order derivative using the skewed right-sided difference, skewed left-sided difference and central difference scheme. And to write a MATLAB program to evaluate the second-order derivative of the analytical function `e^x*cosx` and compare it with the 3 numerical approximations derived earlier.
2. Derivation of the approximation schemes:
- Skewed right-sided approximation scheme:
In skewed right-sided scheme, to compute the value of the derivative at a point, we take into consideration the point where we want to find derivative and points to its right side. For finding derivative of nth order which is accurate upto order 'm' using skewed scheme, number of points required is m+n. Since here we have to find second order derivative which is fourth order accurate, total points used in the scheme is 2+4=6. The points will be fi, fi+1, fi+2, fi+3, fi+4, fi+5. Now we write a general equation for second order derivative by assigning some coefficients to each of these six terms. And then we will try to obtain equations and solve them to find out the value of these coefficients.
General form of the scheme:
`(del^2f)/(delx^2) = a*f(i)+b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+g*f(i+5)` ................................................ Eq. (1)
Now, we can apply Taylor series expansion to each of the six terms on the right side as follows:
`a*f(i) = a*f(i)`
`b*f(i+1)=b*f(i)+b*(f')/(1!)(1dx)+b*(f'')/(2!)(1dx)^2+b*(f''')/(3!)(1dx)^3+b*(f'''')/(4!)(1dx)^4+b*(f''''')/(5!)(1dx)^5+b*(f'''''')/(6!)(1dx)^6+...`
`c*f(i+2)=c*f(i)+c*(f')/(1!)(2dx)+c*(f'')/(2!)(2dx)^2+c*(f''')/(3!)(2dx)^3+c*(f'''')/(4!)(2dx)^4+c*(f''''')/(5!)(2dx)^5+c*(f'''''')/(6!)(2dx)^6+...`
`d*f(i+3)=d*f(i)+d*(f')/(1!)(3dx)+d*(f'')/(2!)(3dx)^2+d*(f''')/(3!)(3dx)^3+d*(f'''')/(4!)(3dx)^4+d*(f''''')/(5!)(3dx)^5+d*(f'''''')/(6!)(3dx)^6+...`
`e*f(i+4)=e*f(i)+e*(f')/(1!)(4dx)+e*(f'')/(2!)(4dx)^2+e*(f''')/(3!)(4dx)^3+e*(f'''')/(4!)(4dx)^4+e*(f''''')/(5!)(4dx)^5+e*(f'''''')/(6!)(4dx)^6+...`
`g*f(i+5)=g*f(i)+g*(f')/(1!)(5dx)+g*(f'')/(2!)(5dx)^2+g*(f''')/(3!)(5dx)^3+g*(f'''')/(4!)(5dx)^4+g*(f''''')/(5!)(5dx)^5+g*(f'''''')/(6!)(5dx)^6+...`
Adding the above written six equations and grouping terms of the same order of derivative in RHS, we get:
`a*f(i)+b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+g*f(i+5)`
` = (a+b+c+d+e+g)*f+(b+2c+3d+4e+5g)*f'dx+(1/2b+4/2c+9/2d+16/2e+25/2g)*f''(dx)^2`
`+(1/6b+8/6c+27/6d+64/6e+125/6g)*f'''(dx)^3+(1/24b+16/24c+81/24d+256/24e+625/24g)*f''''(dx)^4`
`+(1/120b+32/120c+243/120d+1024/120e+3125/120g)*f'''''(dx)^5+(1/720b+64/720c+729/720d+4096/720e+15625/720g)*f''''''(dx)^6+...`
This above Taylor series expansion is then substituted in Equation (1) and the coefficients of the different order derivative of f are compared. On the LHS of equation (1) the coefficient of `f''` is one and coefficients of all other derivatives of f are zero. So we get the following equations after comparing the coefficients:
`a+b+c+d+e+g=0`
`b+2c+3d+4e+5g=0`
`(1/2b+4/2c+9/2d+16/2e+25/2g)*(dx)^2=1`
`b+8c+27d+64e+125g=0`
`b+16c+81d+256e+625g=0`
`b+32c+243d+1024e+3125g=0`
`(b+64c+729d+4096e+15625g)/720*dx^6` ------------------> This last equation only represents the coefficient of the error term. It is not necessarily equal to zero. And it will not be used for solving the above set of equations. The set of six equations above can be represented using matrix multiplication:
`[[1,1,1,1,1,1],[0,1,2,3,4,5],[0,1/2dx^2,4/2dx^2,9/2dx^2,16/2dx^2,25/2dx^2],[0,1,8,27,64,125],[0,1,16,81,256,625],[0,1,32,243,1024,3125]]*[[a],[b],[c],[d],[e],[g]]=[[0],[0],[1],[0],[0],[0]]`
Solving this matrix, we will get the values of a,b,c,d,e,g:
`a=45/(12dx^2)`
`b=-154/(12dx^2)`
`c=214/(12dx^2)`
`d=-156/(12dx^2)`
`e=61/(12dx^2)`
`g=-10/(12dx^2)`
Substituting these values in Equation (1), we get the skewed right-sided approximation scheme as:
`(del^2f)/(delx^2) = (45*f(i)-154*f(i+1)+214*f(i+2)-156*f(i+3)+61*f(i+4)-10*f(i+5))/(12dx^2)`
Now, we have to prove that this scheme is fourth order accurate. In order to do so, we will find have to find the coefficient of the error term `(f'''''')`. Its coefficient is `(b+64c+729d+4096e+15625g)/720*dx^6`. After substituting the values of a,b,c,d,e,g, we get the coefficient as `-548dx^4`. Since the coefficient is non-zero and the power of dx term is 4, so the scheme is fourth order accurate. - Skewed left-sided approximation scheme:
In skewed left-sided scheme, to compute the value of the derivative at a point, we take into consideration the point where we want to find the derivative and points to its left side. For finding derivative of nth order which is accurate upto order 'm' using skewed scheme, number of points required is m+n. Since here we have to find second order derivative which is fourth order accurate, total points used in the scheme is 2+4=6. The points will be fi, fi-1, fi-2, fi-3, fi-4, fi-5. Now we write a general equation for second order derivative by assigning some coefficients to each of these six terms. And then we will try to obtain equations and solve them to find out the value of these coefficients.
General form of the scheme:
`(del^2f)/(delx^2) = a*f(i)+b*f(i-1)+c*f(i-2)+d*f(i-3)+e*f(i-4)+g*f(i-5)` ................................................ Eq. (2)
Now, we can apply Taylor series expansion to each of the six terms on the right side as follows:
`a*f(i) = a*f(i)`
`b*f(i-1)=b*f(i)-b*(f')/(1!)(1dx)+b*(f'')/(2!)(1dx)^2-b*(f''')/(3!)(1dx)^3+b*(f'''')/(4!)(1dx)^4-b*(f''''')/(5!)(1dx)^5+b*(f'''''')/(6!)(1dx)^6+...`
`c*f(i-2)=c*f(i)-c*(f')/(1!)(2dx)+c*(f'')/(2!)(2dx)^2-c*(f''')/(3!)(2dx)^3+c*(f'''')/(4!)(2dx)^4-c*(f''''')/(5!)(2dx)^5+c*(f'''''')/(6!)(2dx)^6+...`
`d*f(i-3)=d*f(i)-d*(f')/(1!)(3dx)+d*(f'')/(2!)(3dx)^2-d*(f''')/(3!)(3dx)^3+d*(f'''')/(4!)(3dx)^4-d*(f''''')/(5!)(3dx)^5+d*(f'''''')/(6!)(3dx)^6+...`
`e*f(i-4)=e*f(i)-e*(f')/(1!)(4dx)+e*(f'')/(2!)(4dx)^2-e*(f''')/(3!)(4dx)^3+e*(f'''')/(4!)(4dx)^4-e*(f''''')/(5!)(4dx)^5+e*(f'''''')/(6!)(4dx)^6+...`
`g*f(i-5)=g*f(i)-g*(f')/(1!)(5dx)+g*(f'')/(2!)(5dx)^2-g*(f''')/(3!)(5dx)^3+g*(f'''')/(4!)(5dx)^4-g*(f''''')/(5!)(5dx)^5+g*(f'''''')/(6!)(5dx)^6+...`
Adding the above written six equations and grouping terms of the same order of derivative in RHS, we get:
`a*f(i)+b*f(i-1)+c*f(i-2)+d*f(i-3)+e*f(i-4)+g*f(i-5)`
` = (a+b+c+d+e+g)*f-(b+2c+3d+4e+5g)*f'dx+(1/2b+4/2c+9/2d+16/2e+25/2g)*f''(dx)^2`
`-(1/6b+8/6c+27/6d+64/6e+125/6g)*f'''(dx)^3+(1/24b+16/24c+81/24d+256/24e+625/24g)*f''''(dx)^4`
`-(1/120b+32/120c+243/120d+1024/120e+3125/120g)*f'''''(dx)^5+(1/720b+64/720c+729/720d+4096/720e+15625/720g)*f''''''(dx)^6+...`
This above Taylor series expansion is then substituted in Equation (2) and the coefficients of the different order derivative of f are compared. On the LHS of equation (2) the coefficient of `f''` is one and coefficients of all other derivatives of f are zero. So we get the following equations after comparing the coefficients:
`a+b+c+d+e+g=0`
`-(b+2c+3d+4e+5g)=0`
`(1/2b+4/2c+9/2d+16/2e+25/2g)*(dx)^2=1`
`-(b+8c+27d+64e+125g)=0`
`b+16c+81d+256e+625g=0`
`-(b+32c+243d+1024e+3125g)=0`
`(b+64c+729d+4096e+15625g)/720*dx^6` ------------------> This last equation only represents the coefficient of the error term. It is not necessarily equal to zero. And it will not be used for solving the above set of equations. The set of six equations above can be represented using matrix multiplication:
`[[1,1,1,1,1,1],[0,1,2,3,4,5],[0,1/2dx^2,4/2dx^2,9/2dx^2,16/2dx^2,25/2dx^2],[0,1,8,27,64,125],[0,1,16,81,256,625],[0,1,32,243,1024,3125]]*[[a],[b],[c],[d],[e],[g]]=[[0],[0],[1],[0],[0],[0]]`
Solving this matrix, we will get the values of a,b,c,d,e,g:
`a=45/(12dx^2)`
`b=-154/(12dx^2)`
`c=214/(12dx^2)`
`d=-156/(12dx^2)`
`e=61/(12dx^2)`
`g=-10/(12dx^2)`
Substituting these values in Equation (2), we get the skewed left-sided approximation scheme as:
`(del^2f)/(delx^2) = (45*f(i)-154*f(i-1)+214*f(i-2)-156*f(i-3)+61*f(i+4)-10*f(i-5))/(12dx^2)`
Now, we have to prove that this scheme is fourth order accurate. In order to do so, we will find have to find the coefficient of the error term `(f'''''')`. Its coefficient is `(b+64c+729d+4096e+15625g)/720*dx^6`. After substituting the values of a,b,c,d,e,g, we get the coefficient as `-548dx^4`. Since the coefficient is non-zero and the power of dx term is 4, so the scheme is fourth order accurate. - Central Difference Scheme:
In central difference scheme, to compute the value of the derivative at a point, we take into consideration the point where we to find derivative and points on both sides of that point. For finding derivative of nth order which is accurate upto order 'm' using central difference scheme, number of points required is m+n-1. Since here we have to find second order derivative which is fourth order accurate, total points used in the scheme is 2+4-1=5. The points will be fi-2, fi-1 , fi, fi+1, fi+2. Now we write a general equation for second order derivative by assigning some coefficients to each of these five terms. And then we will try to obtain equations and solve them to find out the value of these coefficients.
General form of the scheme:
`(del^2f)/(delx^2) = a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2)` ........................................................................... Eq. (3)
Now, we can apply Taylor series expansion to each of the six terms on the right side as follows:
`a*f(i-2) = a*f(i)-a*(f')/(1!)(2dx)+a*(f'')/(2!)(2dx)^2-a*(f''')/(3!)(2dx)^3+a*(f'''')/(4!)(2dx)^4-a*(f''''')/(5!)(2dx)^5+...`
`b*f(i-1)=b*f(i)-b*(f')/(1!)(1dx)+b*(f'')/(2!)(1dx)^2-b*(f''')/(3!)(1dx)^3+b*(f'''')/(4!)(1dx)^4-b*(f''''')/(5!)(1dx)^5+...`
`c*f(i)=c*f(i)`
`d*f(i+1)=d*f(i)+d*(f')/(1!)(1dx)+d*(f'')/(2!)(1dx)^2+d*(f''')/(3!)(1dx)^3+d*(f'''')/(4!)(1dx)^4+d*(f''''')/(5!)(1dx)^5+...`
`e*f(i+2)=e*f(i)+e*(f')/(1!)(2dx)+e*(f'')/(2!)(2dx)^2+e*(f''')/(3!)(2dx)^3+e*(f'''')/(4!)(2dx)^4+e*(f''''')/(5!)(2dx)^5+...`
Adding the above written five equations and grouping terms of the same order of derivative in RHS, we get:
`a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2)`
` = (a+b+c+d+e)*f+(-2a-b+d+2e)*f'dx+(4/2a+1/2b+1/2d+4/2e)*f''(dx)^2`
`+(-8/6a-1/6b+1/6d+8/6e)*f'''(dx)^3+(16/24a+1/24b+1/24d+16/24e)*f''''(dx)^4`
`+(-32/120a-1/120b+1/120d+32/120e)*f'''''(dx)^5+(64/720a+1/720b+1/720d+64/720e)*f''''''(dx)^6+...`
This above Taylor series expansion is then substituted in Equation (3) and the coefficients of the different order derivative of f are compared. On the LHS of equation (3) the coefficient of `f''` is one and coefficients of all other derivatives of f are zero. So we get the following equations after comparing the coefficients:
`a+b+c+d+e=0`
`-2a-b+d+2e=0`
`(4/2a+1/2b+1/2d+4/2e)*(dx)^2=1`
`-8a-b+d+8e=0`
`16a+b+d+16e=0`
`(-32a-b+d+32e)/120*dx^5`
`(64a+b+d+64e)/720*dx^6`
The last two equations only represents the coefficient of the error terms. They are not necessarily equal to zero. And they will not be used for solving the above set of equations. The set of five equations above can be represented using matrix multiplication:
`[[1,1,1,1,1],[-2,-1,0,1,2],[4/2dx^2,1/2dx^2,0,1/2dx^2,4/2dx^2],[-8,-1,0,1,8],[16,1,0,1,16]]*[[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`
Solving this matrix, we will get the values of a,b,c,d,e:
`a=-1/(12dx^2)`
`b=16/(12dx^2)`
`c=-30/(12dx^2)`
`d=16/(12dx^2)`
`e=-1/(12dx^2)`
Substituting these values in Equation (3), we get the central difference approximation scheme as:
`(del^2f)/(delx^2) = (-1*f(i-2)+16*f(i-1)-30*f(i)+16*f(i+1)-1*f(i+2))/(12dx^2)`
Now, we have to prove that this scheme is fourth order accurate. In order to do so, we will find have to find the coefficient of the error term `(f''''')`. Its coefficient is `(-32a-b+d+32e)/120*dx^5`. After substituting the values of a,b,c,d,e, we get the coefficient as `0`. Since the coefficient is zero there is no term of that order, so we check the coefficient of the next error term `(f'''''')`. Its coefficient is `(64a+b+d+64e)/720dx^6`. After substituting the values of a,b,c,d,e, we get the coefficient as `-1/90dx^4`. Since the coefficient is non-zero and the power of dx term is 4, so the scheme is fourth order accurate.
3. Explanation of the MATLAB Code:
- Firstly, we declare the variables x and dx by assigning them values. x contains the value of the point where we want to find the derivative. dx contains the range of values over which the error in the different approximation schemes will be compared. linspace command is used to divide the space between pi/4 and pi/400 into 1000 equal points. The variable x is declared and assigned global scope using keyword global so that it can be accessed inside the functions.
- Next a function has been created in which the function, whose derivative is to be calculated, is defined. It takes a particular value of x as input and return the value of the function at that point.
- A variable called analytical_deriv is declared that contains the analytically found derivative of the function at x.
- We then create three functions, one each for the three different methods: right-skewed, left-skewed and central difference. They are called forward_skewed, backward_skewed and central_diff respectively. These functions accept a value of dx and for that value of dx return the value of the numerical derivative using their respective schemes. The formula for these schemes is the same as derived earlier above.
- A for loop is created which is used to access different values of dx and all three functions defined in the previous steps are called one by one, by passing each value of dx as a parameter in each iteration of the for loop. The returned value from these functions is the numerical derivative. Three row matrices for storing error, one for each type of scheme, are created which contain the absolute difference of the analytical derivative and their respective numerical derivative returned from the function call. These error matrices are named forward_scheme_error, backward_scheme_error and central_scheme_error.
- The errors for all three schemes are plotted using loglog command. A loglog plot is used as the magnitude of error becomes very small for small values of dx and it is difficult to see the variation in a normal plot. The axis labels, title and legend have been mentioned accordingly.
4. MATLAB Code:
clear all
close all
clc
global x=pi/3; %The point where derivative is calculated. Its scope is set to global.
dx=linspace(pi/4,pi/400,1000); %A row matrix containing 1000 equispaced points between pi/4 and pi/400 to represent range of dx(spacing between the discrete points)
%This function returns the function whose derivative is to be calculated
function out = f(u)
out=exp(u)*cos(u);
endfunction
analytical_deriv=-2*exp(x)*sin(x); %Stores analytical derivative of the given function
%This function returns value of right skewed numerical derivative at x.
function fwd=forward_skewed(dx)
global x;
fwd=(45*f(x)-154*f(x+1*dx)+214*f(x+2*dx)-156*f(x+3*dx)+61*f(x+4*dx)-10*f(x+5*dx))/(12*dx^2);
end
%This function returns value of left skewed numerical derivative at x.
function bwd=backward_skewed(dx)
global x;
bwd=(45*f(x)-154*f(x-1*dx)+214*f(x-2*dx)-156*f(x-3*dx)+61*f(x-4*dx)-10*f(x-5*dx))/(12*dx^2);
end
%This function returns value of central difference numerical derivative at x.
function cds=central_diff(dx)
global x;
cds=(-1*f(x-2*dx)+16*f(x-dx)-30*f(x)+16*f(x+dx)-1*f(x+2*dx))/(12*dx^2);
end
%For loop is used to calculate error for each scheme for different values of dx by finding absolute difference of analytical and its respective numerical derivative
for i=1:length(dx)
forward_scheme_error(i)=abs(analytical_deriv-forward_skewed(dx(i)));
backward_scheme_error(i)=abs(analytical_deriv-backward_skewed(dx(i)));
central_scheme_error(i)=abs(analytical_deriv-central_diff(dx(i)));
end
%The error for all three schemes for different values of dx is plotted on a logarithmic scale.
loglog(dx,forward_scheme_error,'r',dx,backward_scheme_error,'g',dx,central_scheme_error,'b');
xlabel('Range of dx'), ylabel('Error Magnitude')
title('Error comparison of three methods of approximation over a range of dx')
legend('Error Forward(Right) Skewed','Error Backward(Left) Skewed','Error Central Difference','location','northwest')
axis([pi/400 pi/4]);5. Output Plot:
6. Results and Conclusion:
From the graph, we can see that central difference scheme is more accurate than skewed scheme for most values of dx. And as expected, we can see that error decreases as the value of dx is decreased. If we further decrease the value of dx, after a point error starts increasing again. This is due to round-off and truncation error. This round-off error out-weighs the accuracy achieved through closer spacing between the points and hence the error starts increasing. So an ideal dx value, depending on the problem, needs to be chosen such that error is minimum.
Also, even though the skewed schemes are less accurate than central difference scheme, they are useful in certain situations. Like at the boundaries of a domain, when there are no points on both sides of the node on the border, central difference scheme cannot be applied. So in that case a skewed scheme can be used. Also in problems where the solution at a particular point depends only on upstream conditions, then in such cases the points downstream have no effect on the result at the current point. So the value is only influenced by the value of that property at only upstream points. Hence a skewed scheme will be more effective and give better results.