Study of Taylor Table Method to derive numerical formula for nth order derivative and using Taylor table method to derive formula and then compare the trend of its numerical errors for Central difference and Skewed difference schemes for 4th order approximation of 2nd order derivative
Taylor table method has proven to be effective and less clumsy method to derive numerical schemes for different order of derivatives.In the challenge,it is given to derive central and skewed schemes of 2nd order derivative,so one by one they will be derived and respective functions will also be crafted along with each of them.
Lets write code of functions of input function and its analytical derivative,so that in future there is ease in navigation through the challenge by simply replacing the input function and its derivative by there respective functions
1)Input Function
`f(x)=exp(x)*cos(x)`
and hence code for the function for the input function is
function [output]=inpfun(x)
output=exp(x)*cos(x);
endfunction2) Analytical 2nd order derivative
`f(x)^(\'\')=-2*exp(x)*sin(x)`
and hence the code for its function is
function [output]=analytical_derivative(x)
output=-2*exp(x)*sin(x);
endfunction
So now first,it has to be understood how a taylor table is crafted and how to decide the number of points about the x value to be taken to calculate a derivative at that x of a particular order i.e as 2nd order derivative of 4th order is required, how many points about x to be taken to derive its formula.
It is of paramount importance to understand that the number of points that will be taken to calculate 2nd order or for that matter any order of derivative has to be decided before hand and it depends on the following factors:-
1.)The order of derivative that is calculated by the formula
2.)the order of the approximation of numerical scheme decided
3.)The difference scheme that is selected i.e central or skewed
Lets first understand what taylors series is:-
`f(x+dx)=f(x)+(dx*f(x)^\')/(1!)+(dx^2*f(x)^(\'\'))/(2!)+(dx^3*f(x)^(\'\'\'))/(3!)........`
In this manner there are more terms with increasing powers of dx and increasing order of derivative.
Now certain nomenclature has to be understood to save some time:-
`f(x)=fi`
`f(x+dx)=f(i+1)`
`f(x-dx)=f(i-1)`
`f(x+2*dx)=f(i+2)`...So in general
`f(x+n*dx)=f(i+n)`
So the same taylor equation becomes
`f(i+1)=f(i)+(dx*f(i)^\')/(1!)+(dx^2*f(i)^(\'\'))/(2!)+(dx^3*f(i)^(\'\'\'))/(3!)........``
In the same way if f(i+n) has to be calculated the above taylor equation can be replaced by
`f(i+n)=f(i)+((n*dx)*f(i)^\')/(1!)+((n*dx)^2*f(i)^(\'\'))/(2!)+((n*dx)^3*f(i)^(\'\'\'))/(3!)........`
So now it can be understood that any point around x can be represented in the form of function and its derivatives at x.
So another thing that have to be understood depending on the order of numerical scheme and order of numerical scheme,number of points that have to be taken are decided:-
For example if 4th order approximation of 2nd order derivative is required,then 6 points will be required for the formula.But why it is decided by the order of the approximation and order of derivative,it will be explained after explaining taylor table.
For the time being lets assume we know the number points required for fourth order approximation of second order derivative
So according to the question formula for 4th order approximation of 2nd order derivative required:-
So it is assumed that
`f(x)^(\'\')=f(i)^(\'\')=a*f(i)+b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+g*f(i+5)` where a,b,c,d,e,g are random constants.Now it is important to find out these constants so as to get a formula of 2nd order derivative.Lets call a,b,c,d,e and g as formula coefficients for ease of refering to them.
Note:-All the points that are taken are to the right side of x just to prove a point. These point that are taken can be any points in the vicinity of the x or i^th point that is in consideration.For the explaining purpose,skewed right-sided scheme is solved here from table method as an example.
So the idea is that taylor series is used to expand,f(i),f(i+1),f(i+2),f(i+3),f(i+4) and f(i+5) in the above expression so that we will get a form of equation something like this after simplyfying:-
`f(i)^(\'\')=J*f(i)+K*f(i)^(\')+L*f(i)^(\'\')+M*f(i)^(\'\'\')+N*f(i)^(\'\'\'\')+Q*f(i)^(v)`....i
where J,K,L,M,N,Q are equations of a,b,c,d,e,g along with terms of dx which are formed due to coefficients of terms in taylor series expansion.Now let call J,K,L,M,N and Q as taylor coefficients for ease of refering to them.
The next step is that since derivative terms on both the sides are present hence a comparison of taylor coefficients of derivative terms can be made between LHS and RHS.
Citing above logic,it can be inferred:
J=K=M=N=Q=0 and L=1
Since second derivative term is the only derivative term present on the LHS and has one as coefficient.
As it can be seen that J=0,K=0,L=1,M=0,N=0 and Q=0 are a set of equations in variables a,b,c,d,e and g, it can get cumbersome to draw out these equations one by one and solve them. Hence Taylor table is formed to ease out the steps and assimilate all the step into one.
The taylor table is as follows:-
This is how a table is formed.To explain it further,from left to right is the respective expansion of taylor series and from top to bottom is the taylor coefficient of the derivatives at i.
From the table, equations become:-
1)`J=a+b+c+d+e+g` and J=0 so by simplyfying
`a+b+c+d+e+g=0`
2)`K=((b+2*c+3*d+4*e+5*g)*dx)/(1!)` and K=0 and since dx can\'t be zero,therefore
`b+2*c+3*d+4*e+5*g=0`
3) `L=((b+4c+9d+16e+25g)*dx^2)/(2!)` and L=1 and by further symplifying we get
`b+4c+9d+16e+25g=2/(dx^2)`
4) `M=((b+8c+27d+64e+125g)*dx^3)/(3!)` and M=0 and since dx can\'t be zero,therefore
`b+8c+27d+64e+125g=0`
5) N=`((b+16c+81d+256e+625g)*dx^4)/(4!)` and N=0 and since dx can\'t be zero,therefore `b+16c+81d+256e+625g=0`
6) `Q=((b+32c+243c+1024e+3125g)*dx^5)/(5!)` and Q=0 and since dx can\'t be zero,therefore `b+32c+243c+1024e+3125g=0`
Solving all the equations in bold will fetch the values of formula coefficients : a,b,c,d.e and g and this will ultimately form a formula for numerical second order derivative.
Now as it is explained what taylor table is and how it is used.
Procedure to determine number of points required to derive formula for a particular order of approximation and order of derivaive
Lets understand why n+2 points are taken for nth order approximation of 2nd order derivative so that if in future we have to find out how many points are required for a numerical scheme we are able to do it.But before that we require two things to know:-
1)order of approximation
2)order of derivative in calculation
Lets take this situation and assume that it is not known how many points are to be taken for 4th order approximation 2nd order derivative
First step:-
We know that after representing 2nd order derivative in terms of a polynomial of points about point x, we are ultimately going to get a equation of derivatives with taylor coefficients like this:-
`f(i)^(\'\')=J*f(i)+K*f(i)^(\')*dx+L*f(i)^(\'\')*dx^2+M*f(i)^(\'\'\')*dx^3+N*f(i)^(\'\'\'\')*dx^4+Q*f(i)^(v)*dx^5+R*f(i)^(vi)*dx^6+S*f(i)^(vii)*dx^7`
.....ii
Its just that right now it is not known till how many terms the above equation has to be written.Now here we have taken till dx^7 and it can be beyond that as well but for the sake of space it is taken till that. Now this equation is similar to equation (i) which was derived earlier knowing the number of points but here number of points is not known.
One more difference is that dx terms are written seperately from coefficients.Previously dx terms were also part of taylor coefficients like K or L.
Second step:- Now since 2nd order derivative is desired, divide every term of equation (ii) by dx^2....(If it is nth order derivative, divide by dx^n ) .....we get
`f(i)^(\'\')=(J*f(i))/(dx^2)+(K*f(i)^(\'))/(dx)+L*f(i)^(\'\')+M*f(i)^(\'\'\')*dx+N*f(i)^(\'\'\'\')*dx^2+Q*f(i)^(v)*dx^3+R*f(i)^(vi)*dx^4+S*f(i)^(vii)*dx^5`
......iii
Third Step:- Now since 4th order approximation is desired,that means the error term must contain dx^4. Therefore dx^4 term in equation (iii) is the error term and the terms after that with dx with higher power than are four are omitted and considered as negligible.(Note:- if nth order approximation is required then term with x^n becomes error term and terms with higher power that n become negligible and are therefore neglected)
So the actual formula contains the terms till dx^3 and the term dx^4 is considered negligible but is still important to indicate the order of numerical approximation.After doing this,the equation we get is
`f(i)^(\'\')=(J*f(i))/(dx^2)+(K*f(i)^(\'))/(dx)+L*f(i)^(\'\')+M*f(i)^(\'\'\')*dx+N*f(i)^(\'\'\'\')*dx^2+Q*f(i)^(v)*dx^3+R*f(i)^(vi)*dx^4`
....iv
The term with taylor coefficient R denotes the error term but it is actually omitted from the actual formula and so actual formula becomes
`f(i)^(\'\')=(J*f(i))/(dx^2)+(K*f(i)^(\'))/(dx)+L*f(i)^(\'\')+M*f(i)^(\'\'\')*dx+N*f(i)^(\'\'\'\')*dx^2+Q*f(i)^(v)*dx^3`
.....v
Fourth Step:- Now count the number of terms in the equation (v). We can observe that in the equation there are 6 terms and now lets rewrite equation (ii) till six terms:-
`f(i)^(\'\')=J*f(i)+K*f(i)^(\')*dx+L*f(i)^(\'\')*dx^2+M*f(i)^(\'\'\')*dx^3+N*f(i)^(\'\'\'\')*dx^4+Q*f(i)^(v)*dx^5`
Fifth Step:- Now if there n terms in the equation obtained in the fourth step, then n points are required in the vicinity of x, including it as well. Now since 6 terms are present in taylor table formula in the fourth step,therefore six points are required in the vicinity of x and around x.
The reason for this is:-
It is established that taylor table equation will have 6 terms.Therefore therefore there will be 6 taylor constants which are namely:-J,K,L,M,N and Q. These constants are basically equations which must have 6 Formula constants defining each one of them in order for them to be solvable.And since there are 6 formula constants (a,b,c,d,e and g),that means there are 6 points in the formula.
Special case:-Some times it happens due to a particular selection of points renders value of one of the variable equal to zero and in this case it may happen to variable g when central difference scheme is selected.This is construed as the same order of approximation can be done with one less point.But it is special case of central difference scheme but other than this all other manner of selection of points will lead to selection of 6 point for fourth order approximation of 2nd order derivative.But no matter how the selection of points is done,follow the above procedure for number of points. How to counter central differencing scheme will be told at a later stage.
Now that the methodology ,to decide the number of points to be selected, is established. Taylor table method can be applied to get the formula for the numerical derivative
Lets derive formula for numerical schemes given in the challenge
1.)Skewed right-sided difference scheme
This scheme involves taking points all on the right side of x including x itself.
lets rewrite
`f(x)^(\'\')=f(i)^(\'\')=a*f(i)+b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+g*f(i+5)`
This scheme,which was taken as example for explanation, was solved above and equations that were obtained are:-
`a+b+c+d+e+g=0`
`b+2*c+3*d+4*e+5*g=0`
`b+4c+9d+16e+25g=2/(dx^2)`
`b+8c+27d+64e+125g=0`
`b+16c+81d+256e+625g=0`
`b+32c+243c+1024e+3125g=0`
To solve these equations,assume value of dx=1 and when all the values of formula variables are obtained,then simply divide each one of them by dx^2
Solving these equation by matrix inversion method we get
The code
clear all
close all
clc
a=[1 1 1 1 1 1;0 1 2 3 4 5;0 1 4 9 16 25;0 1 8 27 64 125;0 1 16 81 256 625;0 1 32 243 1024 3125]
c=[0;0;2;0;0;0]
x=inv(a)*cThe output
The formula variable values were converted from decimals to fraction and divided by dx^2,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)`
Now the last thing that has to be done is is check the error term.As it can be pointed out that there was an error term which had R as its Taylor coefficient.It has to be checked whether R is zero or not.If R is zero,then the approximation will not be of the 4th order.
(Note:-R is considered zero in the formula as a consequense of dx^4 to be negligible but the sum of the coeffcients in the R should not amount to zero for the whole formula to 4th order approximation)
Error term is the one on the sixth coloumn which is quite apparent from the table above
And hence R coeffcient is
`R=b+64c+729d+4096e+15625g`
Value of R is
R=-548...Solved in the matlab
Since Value of R is not zero therefore The formula is of the fourth order.
So Now final formula becomes
`f(i)^(\'\')=(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 its important put this formula into a function and then find out error in this formula by subtracting any value from it by analytical derivative of input function.Here in this formula f(i) will be replaced by inpfun(x).
The code for it function is:-
function [output]=right_side_fourth_order_error(x,dx)
numerical_derivative=(45*inpfun(x)-154*inpfun(x+dx)+214*inpfun(x+2*dx)-156*inpfun(x+3*dx)+61*inpfun(x+4*dx)-10*inpfun(x+5*dx))/(12*dx^2);
output=abs(numerical_derivative-analytical_derivative(x));
endfunction Since skewed right sided difference scheme was taken as example to explain the procedure, steps to solve were not that distinct. But now every scheme,hence forward, would follow a procedure with steps to further clarify how the schemes have to be followed
2)Skewed left sided difference scheme:-
Step 1:- Since this is 4th order approximation of 2nd order derivative, therefore 6 points will considered as explained before.
The six points are:-
i,i-1,i-2,i-3,i-4 and i-5
Step 2:-
The equation:-
`f(i)^(\'\')=a*f(i)+b*f(i-1)+c*f(i-2)+d*f(i-3)+e*f(i-4)+g*f(i-5)`
Step3:-
Now the taylor table will be formed and it will include error coefficient as well to reduce the steps further:-
Step 4:-
Since the formula is for f(i)\'\',therefore L=1 and rest all are zero except R i.e. J=K=M=N=Q=0
Therefore equations are:-
1)`J=a+b+c+d+e+g` and J=0 so by simplyfying
`a+b+c+d+e+g=0`
2)`K=((-b-2*c-3*d-4*e-5*g)*dx)/(1!)` and K=0 and since dx can\'t be zero,therefore
`-b-2*c-3*d-4*e-5*g=0`
3)`L=((b+4c+9d+16e+25g)*dx^2)/(2!)` and L=1 and by further symplifying we get
`b+4c+9d+16e+25g=2/(dx^2)`
4)`M=((-b-8c-27d-64e-125g)*dx^3)/(3!)`and M=0 and since dx can\'t be zero,therefore
`-b-8c-27d-64e-125g=0`
5)`N=((b+16c+81d+256e+625g)*dx^4)/(4!)`and N=0 and since dx can\'t be zero,therefore `b+16c+81d+256e+625g=0`
6)`Q=((-b-32c-243c-1024e-3125g)*dx^5)/(5!)` and Q=0 and since dx can\'t be zero,therefore `-b-32c-243c-1024e-3125g=0`
And the Taylor coefficient of the error term=`R=b+64c+729d+4096e+15625g`
Step 5:-
Solving the above equations using matrix inversion method in matlab:-
The code:-
clear all
close all
clc
a=[1 1 1 1 1 1;0 -1 -2 -3 -4 -5;0 1 4 9 16 25;0 -1 -8 -27 -64 -125;0 1 16 81 256 625;0 -1 -32 -243 -1024 -3125]
c=[0;0;2;0;0;0]
x=inv(a)*c
R=x(2)+64*x(3)+729*x(4)+4096*x(5)+15625*x(6)And the output is
As it can be observed,the Taylor coefficient of error term i.e R is not equal to zero therefore the formula will be of 4th order approximation.
Values of the formula coefficients were converted from decimal to fraction and divided by dx^2 and there denominators were also equalized by l.c.m method and they were calculated to be:-
`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)`
step 6:-
Substitute in the formula
`f(i)^(\'\')=(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)`
This is the final formula
Now a function will be written in order to calculate error generated from numerical approximation from this formula by subtracting its value from the value gooten from analytical derivative.
The code:-
function [output]=left_side_fourth_order_error(x,dx)
numerical_derivative=(45*inpfun(x)-154*inpfun(x-dx)+214*inpfun(x-2*dx)-156*inpfun(x-3*dx)+61*inpfun(x-4*dx)-10*inpfun(x-5*dx))/(12*dx^2);
output=abs(numerical_derivative-analytical_derivative(x));
endfunction
3)Central difference scheme
Central difference scheme is the one which causes some disparity from other difference schemes given the order of approximation and order of derivatives is the same among the schemes.But it will be pointed out later in the procedure.For now steps will be done in the same way they were done for previous schemes.
Step 1:-
Points for central difference scheme are always selected in odd numbers i.e
i-1,i,i+1 - 3 points
i-2,i-1,i,i+1,i+2 - 5 points
Since this is 4th order approximation of 2nd order derivative, therefore it will be assumed that it is not known that 5 points will give 4th order approximation and it will be assumed that since procedure set for number of points required for 4th order approximation for 2nd order derivative is six,therefore central difference will be set in accordance with that
For now selection of points will be in this manner since 6 points are desired:-
i-2,i-1,i,i+1,i+2 and i+3 (this selection is little skewed and not perfectly central difference).If fourth order approximation 2nd order derivative allow central difference scheme,then i+3 point will automatically get deleted because [i-2,i-1,i,i+1,i+2,i+3 ] and [i-2,i-1,i,i+1,i+2] would not tantamount to same formula unless formula coefficient for f(i+3) gets zero.If coefficient to f(i+3) gets zero only then it is central difference scheme or otherwise central difference scheme is not possible for fourth order approximation 2nd order derivative.
Step 2:-
The equation:-
`f(i)^(\'\')=a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2)+g*f(i+3)`
Step3:-
Now the taylor table will be formed and it will include error coefficient as well to reduce the steps further:-
Step 4:-
Since the formula is for f(i)\'\',therefore L=1 and rest all are zero except R i.e. J=K=M=N=Q=0
Therefore equations are:-
1)`J=a+b+c+d+e+g` and J=0 so by simplyfying
`a+b+c+d+e+g=0`
2)`K=((-2a-b+d+2e+3g)*dx)/(1!)` and K=0 and since dx can\'t be zero,therefore
`-2a-b+d+2e+3g=0`
3)`L=((4a+b+d+4e+9g)*dx^2)/(2!)` and L=1 and by further symplifying we get
`4a+b+d+4e+9g=2/(dx^2)`
4)`M=((-8a-b+d+8e+27g)*dx^3)/(3!)` and M=0 and since dx can\'t be zero,therefore
`-8a-b+d+8e+27g=0`
5)`N=((16a+b+d+16e+81g)*dx^4)/(4!)` and N=0 and since dx can\'t be zero,therefore `16a+b+d+16e+81g=0`
6)`Q=((-32a-b+d+32e+243g)*dx^5)/(5!)` and Q=0 and since dx can\'t be zero,therefore `-32a-b+d+32e+243g=0`
And the coefficient of the error term=`R=64a+b+d+64e+729g`
Step 5:-
Solving the above equations using matrix inversion method in matlab:-
The code:-
clear all
close all
clc
a=[1 1 1 1 1 1;-2 -1 0 1 2 3;4 1 0 1 4 9;-8 -1 0 1 8 27;16 1 0 1 16 81;-32 -1 0 1 32 243]
c=[0;0;2;0;0;0]
x=inv(a)*c
R=64*x(1)+x(2)+x(4)+64*x(5)+729*x(6)And the output is
As it can be observed,the coefficient of error term i.e R is not equal to zero therefore the formula will be of 4th order approximation.
Values of the variables were converted from decimal to fraction and divided by dx^2 and there denominators were also equalized by l.c.m method and they were calculated to be:-
`a=(-1)/(12dx^2)`
`b=(16)/(12*dx^2)`
`c=(-30)/(12*dx^2)`
`d=(16)/(12*dx^2)`
`e=(-1)/(12dx^2)`
`g=-3.0840*10^-18=0`
Step 6:-
We will get the final formula by substituting values for the variables in the equation we started in second step
`f(i)^(\'\')=(-1*f(i-2)+16*f(i-1)-30*f(i)+16*f(i+1)-1*f(i+2))/(12*dx^2)`
There are two observations that can be made from the formula:-
1) the formula coefficient of f(i+3) or g came as ZERO
2) the formula coefficients are symmetric about the i point....i.e
a.) coefficient of f(i-2)=coefficient of f(i+2)
a=e
b.) coefficient of f(i-1)=coefficient of f(i+1)
b=d
So an important conclusion can be made here that
Now since as g or coefficient for f(i+3) has come out as zero,therefore the above derived formula can be called central difference 4th order approximation of 2nd order derivative. But this means the same formula could have been derived with 5 points instead of six.If five points would have been taken then taylor table would have looked like this:-
When only five points were taken the row with g*f(i+3) was deleted. So there are five points and hence five columns are needed to make the equations obtained from them are solvable.Therefore column till f(i)\'\'\'\' was taken.Now naturally the column of `f(i)^(v)` would have become as error term and Taylor coefficient at the fifth column was Q as can be seen in previous taylor tables.Moreover order of approximation would have become as three rather than 4. But since formula for Q without g(as sixth point is not there ):-
`Q=((-32a-b+d+32e)*dx^5)/(5!)` and as previously observed a=e and b=d,
Q=0(not as a consequence of dx^5 or comparison between RHS and LHS,But as a consequence of the formula symmetricity existing in the formula of Q due to a=e and b=d) and hence if error term coeffcient comes as zero ,the column cannot be that of the error term
Therefore this column with `f(i)^(v)` was deleted as Q=0
Therefore column of error term became one column after that i.e column of f(i)^(vi) and Taylor coefficient of error term became R same as every numerical scheme solved before this and R was found to be non zero. Therefore this makes R term Taylor coefficient of error term for this numerical scheme and order of approximation remained 4th order with 5 points instead of 6
Important Note:-
It can be found that as Q automatically came to zero, that means with central differencing one cannot get 3rd order approximation or for that matter any odd numered order of approximation for 2nd order of derivative.
Important tip:-
While trying to find number of points required to derive formula for n order approximation of m order derivative using central differencing scheme, proceed with the procedure to find points as told previously.If the number of points come as odd,then there is no problem in selection of points in the central differencing form but if it comes as even then following steps have to be followed:-
Step 1:So number of points required is an even number equal to n
Step 2:That means n-1 is odd number and therefore select points in central differencing form as it would be done for an odd number.
If n-1 =3
points i-1,i,i+1
n-1=5
points i-2,i-1,i,i+1,i+2
Step 3:Now select one point more than the points which were selected in the previous step
If n-1 =3 and therefore n=4
points i-1,i,i+1,i+2
If n-1=5 and therefore n=6
points i-2,i-1,i,i+1,i+2,i+3
Now once above procedure is carried out, apply taylor table method to the selected points.Now if formula coefficient for the extra point(f(i+3) for n=6 and f(i+2) for n=4) comes as zero after entire variable calculation, only then the formula can be assigned as central differencing otherwise it will remain slightly skewed towards right.After it has been assigned as central differencing of mth approximation then the Taylor coefficient of error term has to be checked whether it amounts to zero or not to verify the order of approximation.This is the procedure that was applied for deriving formula for central differencing scheme for 4th order numerical approximation of 2nd order.
Now as formula for central differencing scheme is derived, the function for its error with the analytical term is
function [output]=central_differ_fourth_order_error(x,dx)
numerical_derivative=(-1*inpfun(x-2*dx)+16*inpfun(x-dx)-30*inpfun(x)+16*inpfun(x+dx)-inpfun(x+2*dx))/(12*dx^2);
output=abs(numerical_derivative-analytical_derivative(x));
endfunction
Comparison of the error from central,skewed left and right numerical schemes for 4th order approximation for 2nd order derivative
Now that the formulas of different numerical schemes have been derived and code for functions for the error generated by them has been written, a methodology have to be devised so as to plot a graph comparing these errors against different values of dx.
Graph:-
Since pattern of error has to be observed at low values of dx hence logarithmic plot has to be plotted instead of linear plot.
Range of dx:-
For the pattern to be trustworthy,it is imperative that input i.e values on the x-axis are uniformly spaced.
Using linspace command to generate points between pi/4000 to pi/4 generates uniformly spaced points but for the linear plot and not for logarithmic plot.
It can be conveyed by a logarithmic scatter plot depicting a y=1 graph
From the logarithmic scatter plot it is quite evident that that linspace command used to generate points between pi/4000 to pi/4 will not generate uniformly spaced input points for logarithmic plot.
Hence another procedure will be applied in which dx will be in terms of power of 10 and its powers will be uniformly spaced.That means following code procedure will be applied:-
p=-3:0.1:0
dx=10.^pIn this way it will generate uniformly spaced input to the logarithmic plot.It is conveyed by y=1 graph plotted in logarithmic plot with above said spacing of dx:-
It is quite evident from the above scatter plot that the above procedure of drafting range for dx provides uniform spaced input in logarithmic plot.
Hence it will be used for the study of comparison carried out further.
Methodology:-
A for loop will be generated in which each loop will be associated with a value of dx and in this loop numerical errors,from the function that were coded before, will be calculated for the dx value dedicated for that loop and stored in an array.These arrays will later be plotted in the form of logarithmic plot against dx.
Above methodology in the form of code is as follows:-
clear all
close all
clc
x=pi/3;
p=-3:0.1:0;
dx=10.^p;
for i=1:length(dx)
central_differ_error(i)=central_differ_fourth_order_error(x,dx(i));
right_differ_error(i)=right_side_fourth_order_error(x,dx(i));
left_differ_error(i)=left_side_fourth_order_error(x,dx(i));
endfor
loglog(dx,central_differ_error,\'color\',\'r\')
hold on
loglog(dx,right_differ_error,\'color\',\'b\')
hold on
loglog(dx,left_differ_error,\'color\',\'g\')
xlabel(\'dx\')
ylabel(\'error\')
legend(\'central differ error\',\'right differ error\',\'left differ error\')And the plot from the above code is
Above plot displays the pattern of numerical errors against dx values.
Analysis of the plot:-
This plot is divided into three parts; the instability in dx values near 10^-3, line plot between 10^-2 and 10^-1 and instability of the plot near 10^0.
Lets focus on the line plot of the graph:-
For the line plot part of the graph it can be stated in terms a linear algebraic equation:-
`Y=Slope.X+C`....which is standard equation for a line in coordinate plane
But here Y=log(f(dx)),f(dx)=error function,X=log(dx) and C is some constant and log is logarithm to the base 10
lets substitute these terms in line equation:-
`log(f(dx))=Slope*log(dx)+logK`
.......where C is depicted as logarithm of constant K
`log(f(dx))=log(dx^(slope))+logk`
`log(f(dx))=log(K.dx^(slope))`
.....removing log from both sides
`f(dx)=K.dx^(Slope)`........a)
This is the expression obtained from the logarithmic plot
Lets analyse error function obtained from Taylor series:-
Numerical derivative=Analytical Derivative +error Term
...This is the general expression
Numerical derivative-Analytical derivative=error Term
As R is the coefficient of the error term in previous derivation of numerical schemes,lets assume it to coefficient of error terms in general
Lets assume further numerical scheme is of the n order approximation. Therefore power of dx in error term will be n
Error term=`f(dx)=R.dx^n + O`.....b)
.......here R not only represents sum of the formula variables like a,b,c etc but also derivative term of taylor series that occurs in the error term and O represents terms containing dx of the higher power than n and higher order derivatives but O becomes negligible as dx values are low where numerical schemes are generally operated. After neglecting O,error term becomes:-
`f(dx)=R.dx^n`.....c)
Now if equation a and equation c are compared, it is found that the equation structure is same and there is just change in variables.One equation defines behaviour of equation in graph plot and other equation regarding the same error function is derived from the analytical taylor series equation.
If both equations studied and analysed, it is found that since both equation represent formula for error term,therefore there variables on the RHS can be equalized. In essence
K=R and Slope=n
i.) K=R
From both the equations; a and c respectively,it understood that K is equal to R. But it is found in calculations of formula for derivatives using taylor table method, that R,as Taylor coefficient of error term always turns up negative in number. And since K=R and therefore either logK is not possible or it represents value in complex plane. From the graph,it can be forund that logK in the equation `log(f(dx))=Slope*log(dx)+logK` is Y-intercept for the logarithmic graph of log(f(dx)) and since there is no origin in an logarithmic graph,that means any line graph in an logarithmic graph will intersect Y-axis in a complex plane i.e Y-intercept exists in complex plane. This explanation conforms with theoritical relation K=R.
ii.)Slope=n
By this relation it would mean that slope of the line graph ,obtained by plotting behaviour of derivative formula due to nth order of approximation against dx, will be equal to order of approximation of the formula itself.
Lets verify this
Slope of the graphs at the middle and dominant part of the graph are calculated by this:-
In the above picture the formula for the slopes and the values are given.
Comparison:-
Now since skewed left,skewed right and central difference schemes that were derived are fourth order approximation therefore
n=4
And it is found that
Slope_central_differ_error = 3.9685
Slope_right_differ_error = 4.0300
Slope_left_differ_error = 3.9743
All the three slopes are approximately equal to 4
This confirms the relation obtained between slope and order of approximation as from the above comprison it is found that `slope~~n`
Important Note:-
This means that if ever validity of a formula for numerical derivative has to be verified,the a logarithmic plot of its numerical error and dx must be drafted and then the error of its line graph must be calculated.The plot being a line graph for most part of the plot is also a necessity. If the slope is approximately equal to the order of approximation for which formula was derived,then the formula is genuine.
Instability near 10^0:
As dx value increases and becomes comparable with x value, the O term in the equation b starts to gather momentum:-
`f(dx)=R.dx^n + O`
The O term no longer remains negligible and starts increasing in value compared to to R.dx^n term. And this makes the curve jagged at the end because of the nonlinearity in an logarithmic plot.
Instability near 10^-3
Now since O term is negligible as dx is kept at a very low value near 10^-3 but then also curve becomes unstable and unpredictable near the end.This is because as error starts to decrease near 10^-3,the round off error due to computation by the processor becomes prominant.If we assume computer to be a perfect contraption to carry out calculations, then the error vs dx slope would have been a straight line towards lower values of dx no matter how small dx gets in terms of its value. But processor of a computer is not perfect,it can be accurate only for so long.Therefore because of these round off error, the straight line graph starts to become unpredictable and jagged at 10^-3.It has to be noted that this error is not due to the error relation but due to error in computation by processor.
Comparison of Central difference scheme and skewed schemes:-
From the graph, it can be found that Central difference scheme is much more accurate than skewed schemes:-left and right.Though all the three are of same order of approximation but central difference scheme is the only one with more accuracy. And one more advantage of central difference scheme is that it requires one point less than its other counter part for the same order of approximation. If it is played around a little with central difference scheme for combination of various selection of points and order of derivatives, it is found that order of approximation will always be an even number.So if anywhere a challenge appears to derive formula for odd order of approximation of any order of derivative using central difference scheme, it can be right away stated that an odd order approximation using central differencing is not possible.
Its application:
These numerical schemes are basically applied in finite difference method to solve different types of problems where a physical domain is discretized into various nodes. Now when central differencing scheme will be applied to a discretized domain, then it will so happen that always a finite number of nodes will be left for which this numerical schemes formula could not be applied.If it is considered the formula which was derived for central differencing namely:
`f(i)^(\'\')=(-1*f(i-2)+16*f(i-1)-30*f(i)+16*f(i+1)-1*f(i+2))/(12*dx^2)`
So in accordance to this formula,for the first node there will not be any nodes before the first node.Therefore this formula cannot be applied to first node.Similarly formula cannot be applied to second,third and the last three nodes of the domain. Now it is sometimes suggested that for the nodes this formula could not calculate the derivative,formula from the skewed schemes be used.But as it can be seen the graph, there is a considerable difference in the accuracy of skewed and central which may give rise to error and this error on later stage may aggrandize to a big value and make solution unstable.Therefore use of central difference scheme,inspite of being more accurate,is not recommended.
Now It is true that left skewed scheme will leave few nodes at the start uncalculated and right skewed scheme may leave few nodes at the end uncalculated.But as it can be observed from the graph their accuracy is near equal to each other. Therefore for whatever nodes left skewed scheme cannot calculate derivative,formula of right skewed scheme can be used and vice-versa provided order of approximation for both schemes is same.This will ensure derivative is calculated for each and every node.
Caution:-
Whenever a numerical scheme is to be used in scientific procedure, hence it would be prudent to first draw out its error vs dx logarithmic graph, and find out for what dx interval the logarithmic plot remains a staright line. And use dx value from that interval for the numerical derivative formula for the scientific procedure.