Week 3 - Taylor table method and Matlab code
Our target is to find Finite Difference formulas for 2nd order derivative with 4th order approximations.
A generalized method for p th order derivative with qth order approximations
By using Taylors series, we can get Finite difference formulae for any order derivative with any order of approximations.
Let us see Some generalized way.
`(d^pf)/dx^p=(af_i+bf_(i+1)+cf_(i+2)+df_(i+3)+ef_(i+4)+.............zf_(i+n-1))/h^p+O(h^q)`
where p is Order of Derivative
q is Order of approximation
n is the number of Grid points
h is the step size
a,b,c,d,e,f,...etc are constants
Following conditions always satisfy for the above equation.
- one-sided(either forward or backward) finite approximations n=p+q.
- Both-sided(central scheme)finite approximations n=p+q-1.
To chuck above conditions let us see some basic examples
- 1st order derivative with the first-order approximation(Forward scheme)
`(df)/dx=(f(x+h)-f(x))/h+O(h)`
Now say Order of Derivative p =1
Order of approximation q=1
number of Grid points n=p+q=2 (Those are `f(x),f(x+h)`)
- 1st order derivative with the secound-order approximation(central scheme)
`(df)/dx=(f(x+h)-f(x-h))/(2h)+O(h^2)`
Order of Derivative p =1
Order of approximation q=2
number of Grid points n=p+q-1=2(Those are `f(x+h)`,`f(x-h)`)
Now let's move on to given question
Here they are asking for 4th order approximation of 2nd derivative.
by formula
p=2,q=2 and n value depends on the type of problem like one-sided or central.
- Skewed right-sided difference
2. Skewed left-sided difference
3.Central difference
1.Skewed right-sided difference(Forward Difference):
Order of Derivative p =2
Order of approximation q=4
Number of grid points =p+q=6
Skewed right-sided means we have to take the grid points positive side of i.
here we have to form FDM by using 6 grid points.those are `i,i+1,i+2,i+3,i+4,i+5`
`(d^2f)/dx^2=(af_i+bf_(i+1)+cf_(i+2)+df_(i+3)+ef_(i+4)+ff_(i+5))/h^2+O(h^4)`
multiply both sides with `h^2`
`h^2(d^2f)/dx^2+O(h^6)=(af_i+bf_(i+1)+cf_(i+2)+df_(i+3)+ef_(i+4)+ff_(i+5))`---------------------------->Equation`
Here we have to convert all `f(i+1),f(i+2),f(i+3),f(i+4) and f(i+5)` into Taylers series about `f_(i)`
`bf_(i+1)=b(f_(i)+(h^1)/(1!)((df)/dx)_i+(h^2)/(2!)((d^2f)/dx^2)_i+(h^3)/(3!)((d^3f)/dx^3)_i+(h^4)/(4!)((d^4f)/dx^4)_i+(h^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
`cf_(i+2)=c(f_(i)+((2h)^1)/(1!)((df)/dx)_i+((2h)^2)/(2!)((d^2f)/dx^2)_i+((2h)^3)/(3!)((d^3f)/dx^3)_i+((2h)^4)/(4!)((d^4f)/dx^4)_i+((2h)^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
`d*f_(i+3)=d(f_(i)+((3h)^1)/(1!)((df)/dx)_i+((3h)^2)/(2!)((d^2f)/dx^2)_i+((3h)^3)/(3!)((d^3f)/dx^3)_i+((3h)^4)/(4!)((d^4f)/dx^4)_i+((3h)^5)/(5!)((d^5f)/dx^5)_i+O(h^6)`
`ef_(i+4)=e(f_(i)+((4h)^1)/(1!)((df)/dx)_i+((4h)^2)/(2!)((d^2f)/dx^2)_i+((4h)^3)/(3!)((d^3f)/dx^3)_i+((4h)^4)/(4!)((d^4f)/dx^4)_i+((4h)^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
`f*f_(i+5)=f(f_(i)+((5h)^1)/(1!)((df)/dx)_i+((5h)^2)/(2!)((d^2f)/dx^2)_i+((5h)^3)/(3!)((d^3f)/dx^3)_i+((5h)^4)/(4!)((d^4f)/dx^4)_i+((5h)^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
Now substitute those all formulas in equation A and following results will get.
LHS
`h^2(d^2f)/dx^2+O(h^6)=`
RHS
`((a+b+c+d+e+f)f_(i)+(b+2c+3d+4e+5f)((h)^1)/(1!)((df)/dx)_i+(b+4c+9d+16e+25f)((h)^2)/(2!)((d^2f)/dx^2)_i+(b+8c+27d+64e+125f)((h)^3)/(3!)((d^3f)/dx^3)_i+(b+16c+81d+256e+625f)((h)^4)/(4!)((d^4f)/dx^4)_i+(b+32c+243d+1024e+3125f)((h)^5)/(5!)((d^5f)/dx^5)_i+(.....)O(h^6))`
now compare coefficients on both sides
`(a+b+c+d+e+f)=0`-------------------------------->`(a+b+c+d+e+f)=0`(a)
`(b+2c+3d+4e+5f)((h)^1)/(1!)=0`---------------->`(b+2c+3d+4e+5f)=0` (b)
`(b+4c+9d+16e+25f)((h)^2)/(2!)=h^2`---------->`(b+4c+9d+16e+25f)=2 (c)`
`(b+8c+27d+64e+125f)((h)^3)/(3!)=0`---------->`(b+8c+27d+64e+125f)=0` (d)
`(b+16c+81d+256e+625f)((h)^4)/(4!)=0`------->`(b+16c+81d+256e+625f)=0`(e)
`(b+32c+243d+1024e+3125f)((h)^5)/(5!)=0`--->`(b+32c+243d+1024e+3125f)=0` (f)
`[[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]][[a],[b],[c],[d],[e],[f]]=[[0],[0],[2],[0],[0],[0]]`
`[[a],[b],[c],[d],[e],[f]]=[[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]]^-1[[0],[0],[2],[0],[0],[0]]`
Matlab Code for inverse
***The sym code used for Fractional output.
%x=A^(-1)*b
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]]
b=[[0];[0];[2];[0];[0];[0]];
x=(inv(A)*b)
x=sym(inv(A)*b)%sysm for fractioal values
output values
`[[a],[b],[c],[d],[e],[f]]=[[15/4],[-77/6],[107/6],[-13],[61/12],[-5/6]]`
Now substitute those values in to `(d^2f)/dx^2=(af_i+bf_(i+1)+cf_(i+2)+df_(i+3)+ef_(i+4)+ff_(i+5))/h^2+O(h^4)`
then `(d^2f)/dx^2=(45f_i-154f_(i+1)+214f_(i+2)-156f_(i+3)+61f_(i+4)-10f_(i+5))/(12h^2)`
2.Skewed left-sided difference(Backward Difference):
Order of Derivative p =2
Order of approximation q=4
Number of grid points =p+q=6
Skewed left-sided means we have to take the grid points negative side of i.
here we have to form FDM by using 6 grid points.Those are `i,i-1,i-2,i-3,i-4,i-5`
`(d^2f)/dx^2=(af_i+bf_(i-1)+cf_(i-2)+df_(i-3)+ef_(i-4)+ff_(i-5))/h^2+O(h^4)`
multiply both sides with `h^2`
`h^2(d^2f)/dx^2+O(h^6)=(af_i+bf_(i-1)+cf_(i-2)+df_(i-3)+ef_(i-4)+ff_(i-5))/h^2`
Here we have to convert all`f(i-1),f(i-2),f(i-3),f(i-4) and f(i-5)` into Taylers series about `f_(i)`
`bf_(i-1)=b(f_(i)-(h^1)/(1!)((df)/dx)_i+(h^2)/(2!)((d^2f)/dx^2)_i-(h^3)/(3!)((d^3f)/dx^3)_i+(h^4)/(4!)((d^4f)/dx^4)_i-(h^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
`cf_(i-2)=c(f_(i)-((2h)^1)/(1!)((df)/dx)_i+((2h)^2)/(2!)((d^2f)/dx^2)_i-((2h)^3)/(3!)((d^3f)/dx^3)_i+((2h)^4)/(4!)((d^4f)/dx^4)_i-((2h)^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
`d*f_(i-3)=d(f_(i)-((3h)^1)/(1!)((df)/dx)_i+((3h)^2)/(2!)((d^2f)/dx^2)_i((3h)^3)/(3!)((d^3f)/dx^3)_i+((3h)^4)/(4!)((d^4f)/dx^4)_i-((3h)^5)/(5!)((d^5f)/dx^5)_i+O(h^6)`
`ef_(i-4)=e(f_(i)-((4h)^1)/(1!)((df)/dx)_i+((4h)^2)/(2!)((d^2f)/dx^2)_i-((4h)^3)/(3!)((d^3f)/dx^3)_i+((4h)^4)/(4!)((d^4f)/dx^4)_i-((4h)^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
`f*f_(i-5)=f(f_(i)-((5h)^1)/(1!)((df)/dx)_i+((5h)^2)/(2!)((d^2f)/dx^2)_i-((5h)^3)/(3!)((d^3f)/dx^3)_i+((5h)^4)/(4!)((d^4f)/dx^4)_i-((5h)^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
Now substitute those all formulas in equation A and following results will get.
LHS
`h^2(d^2f)/dx^2+O(h^6)=`
RHS
`((a+b+c+d+e+f)f_(i)+(-b-2c-3d-4e-5f)((h)^1)/(1!)((df)/dx)_i+(b+4c+9d+16e+25f)((h)^2)/(2!)((d^2f)/dx^2)_i+(-b-8c-27d-64e-125f)((h)^3)/(3!)((d^3f)/dx^3)_i+(b+16c+81d+256e+625f)((h)^4)/(4!)((d^4f)/dx^4)_i+(-b-32c-243d-1024e-3125f)((h)^5)/(5!)((d^5f)/dx^5)_i+(.....)O(h^6))`
now compare coefficients on both sides
`(a+b+c+d+e+f)=0` -------------------------------->`(a+b+c+d+e+f)=0`(a)
`(-b-2c-3d-4e-5f)((h)^1)/(1!)=0`---------------->`(b+2c+3d+4e+5f)=0` (b)
`(b+4c+9d+16e+25f)((h)^2)/(2!)=h^2`---------->`(b+4c+9d+16e+25f)=2 (c)
`(-b-8c-27d-64e-125f)((h)^3)/(3!)=0`---------->`(b+8c+27d+64e+125f)=0` (d)
`(b+16c+81d+256e+625f)((h)^4)/(4!)=0`------->`(b+16c+81d+256e+625f)=0`(d)
`(-b-32c-243d-1024e-3125f)((h)^5)/(5!)=0`--->`(b+32c+243d+1024e+3125f)=0` (f)
`[[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]][[a],[b],[c],[d],[e],[f]]=[[0],[0],[2],[0],[0],[0]]`
Matlab Code for inverse
***The sym code used for Fractional output.
%x=A^(-1)*b
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]]
b=[[0];[0];[2];[0];[0];[0]];
x=(inv(A)*b)
x=sym(inv(A)*b)%sysm for fractioal values
output values
`[[a],[b],[c],[d],[e],[f]]=[[15/4],[-77/6],[107/6],[-13],[61/12],[-5/6]]`
Now substitute those values in to
`(d^2f)/dx^2=(af_i+bf_(i-1)+cf_(i-2)+df_(i-3)+ef_(i-4)+ff_(i-5))/h^2+O(h^4)`
then `(d^2f)/dx^2=(45f_i-154f_(i-1)+214f_(i-2)-156f_(i-3)+61f_(i-4)-10f_(i-5))/(12h^2)`
3.Central difference
Order of Derivative p =2
Order of approximation q=4
Number of grid points =p+q-1=5
the Total five grid points are equally distributed about i
so those are `i,i-1,i-2,i+1,i+2,`
`(d^2f)/dx^2=(af_i+bf_(i+1)+cf_(i+2)+df_(i-1)+ef_(i-2))/h^2+O(h^4)`
multiply both sides with `h^2`
`(d^2f)/dx^2+O(h^6)=(af_i+bf_(i+1)+cf_(i+2)+df_(i-1)+ef_(i-2))`
Here we have to convert all`f(i-1),f(i-2),f(i-3),f(i-4) ` into Taylers series about `f_(i)`
`bf_(i+1)=b(f_(i)+(h^1)/(1!)((df)/dx)_i+(h^2)/(2!)((d^2f)/dx^2)_i+(h^3)/(3!)((d^3f)/dx^3)_i+(h^4)/(4!)((d^4f)/dx^4)_i+(h^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
`cf_(i+2)=c(f_(i)+((2h)^1)/(1!)((df)/dx)_i+((2h)^2)/(2!)((d^2f)/dx^2)_i+((2h)^3)/(3!)((d^3f)/dx^3)_i+((2h)^4)/(4!)((d^4f)/dx^4)_i+((2h)^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
`df_(i-1)=d(f_(i)-(h^1)/(1!)((df)/dx)_i+(h^2)/(2!)((d^2f)/dx^2)_i-(h^3)/(3!)((d^3f)/dx^3)_i+(h^4)/(4!)((d^4f)/dx^4)_i-(h^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
`ef_(i-2)=e(f_(i)-((2h)^1)/(1!)((df)/dx)_i+((2h)^2)/(2!)((d^2f)/dx^2)_i-((2h)^3)/(3!)((d^3f)/dx^3)_i+((2h)^4)/(4!)((d^4f)/dx^4)_i-((2h)^5)/(5!)((d^5f)/dx^5)_i+O(h^6))`
Now substitute those all formulas in equation A and following results will get.
LHS
`h^2(d^2f)/dx^2+O(h^6)=`
RHS
`((a+b+c+d+e)f_(i)+(b+2c-1d-2e)((h)^1)/(1!)((df)/dx)_i+(b+4c+d+4e)((h)^2)/(2!)((d^2f)/dx^2)_i+(b+8c-1d-8e)((h)^3)/(3!)((d^3f)/dx^3)_i+(b+16c+1d+16e)((h)^4)/(4!)((d^4f)/dx^4)_i+(b+32c-1d-32e)((h)^5)/(5!)((d^5f)/dx^5)_i+(.....)O(h^6))`
now compare coefficients on both sides
5 unknows and take 5 equations
`(a+b+c+d+e)=0` -------------------------------->`(a+b+c+d+e)=0`(a)
`(b+2c-1d-2e)((h)^1)/(1!)=0`---------------->`(b+2c-1d-2e)=0` (b)
`(b+4c+d+4e)((h)^2)/(2!)=h^2`---------->`(b+4c+d+4e)=2 (c)
`(b+8c-1d-8e)((h)^3)/(3!)=0`---------->`(b+8c-1d-8e)=0` (d)
`(b+16c+1d+16e)((h)^4)/(4!)=0`------->`(b+16c+1d+16e)=0`(e)
`[[1,1,1,1,1],[0,1,2,-1,-2],[0,1,4,+1,+4],[0,1,8,-1,-8],[0,1,16,+1,+16]][[a],[b],[c],[d],[e]]=[[0],[0],[2],[0],[0]]`
`[[a],[b],[c],[d],[e]]=[[1,1,1,1,1],[0,1,2,-1,-2],[0,1,4,+1,+4],[0,1,8,-1,-8],[0,1,16,+1,+16]]^-1[[0],[0],[2],[0],[0]]`
Matlab Code for inverse
***The sym code used for Fractional output.
%x=A^(-1)*b
A=[[1,1,1,1,1];[0,1,2,-1,-2];[0,1,4,+1,+4];[0,1,8,-1,-8];[0,1,16,+1,+16]];
b=[[0];[0];[2];[0];[0]];
x=(inv(A)*b)
x=sym(inv(A)*b)%sysm for fractioal values`[[a],[b],[c],[d],[e]]=[[-5/2],[4/3],[-1/12],[4/3],[-1/12]]`
Now substitute those values in to
`(d^2f)/dx^2=(af_i+bf_(i+1)+cf_(i+2)+df_(i-1)+ef_(i-2))/h^2+O(h^4)`
then `(d^2f)/dx^2=(-30f_i+16f_(i+1)-1f_(i+2)+16f_(i-1)-1f_(i-2))/(12h^2)`
Final formulas
|
Skewed right-sided difference |
`(d^2f)/dx^2=(45f_i-154f_(i+1)+214f_(i+2)-156f_(i+3)+61f_(i+4)-10f_(i+5))/(12h^2)` |
| Skewed left-sided difference | `(d^2f)/dx^2=(45f_i-154f_(i-1)+214f_(i-2)-156f_(i-3)+61f_(i-4)-10f_(i-5))/(12h^2)` |
|
Central difference |
`(d^2f)/dx^2=(-30f_i+16f_(i+1)-1f_(i+2)+16f_(i-1)-1f_(i-2))/(12h^2)` |
MatLab code
Given
Function: `e^x*cos(x)`
`(d^2f)/dx^2=-2*e^x*sin(x)`
clc
clear all
close all
x=pi/4;
dx=linspace(pi/20, pi/4000, 50);%for 50dx values
Double_derivative_of_f=(-2*exp(x)*sin(x));
for i=1:length(dx)
Fouth_order_Right_skewed(i)=((45*f(x)-154*f(x+dx(i))+214*f(x+2*dx(i))-156*f(x+3*dx(i))+61*f(x+4*dx(i))-10*f(x+5*dx(i))))/(12*dx(i)^2);%Fourth_order_Right_Skewed_method
Fourth_order_Left_skewed(i)=((45*f(x)-154*f(x-dx(i))+214*f(x-2*dx(i))-156*f(x-3*dx(i))+61*f(x-4*dx(i))-10*f(x-5*dx(i))))/(12*dx(i)^2);%Fourth_order_Left_Skewed_method
Fouth_order_central(i)=(-30*f(x)+16*f(x+dx(i))-1*f(x+2*dx(i))+16*f(x-dx(i))-1*f(x-2*dx(i)))/(12*dx(i)^2);%Fourth_order _central
%error functions corresponding method
error_Fouth_order_Right_skewed(i)=abs(Double_derivative_of_f-Fouth_order_Right_skewed(i));
error_Fourth_order_Left_skewed(i)=abs(Double_derivative_of_f-Fourth_order_Left_skewed(i));
error_Fouth_order_central(i)=abs(Double_derivative_of_f- Fouth_order_central(i));
end
figure(1)
title('Error variation with order of approx')
semilogy(dx,error_Fouth_order_Right_skewed ,'g')
hold on
semilogy(dx,error_Fourth_order_Left_skewed,'r')
hold on
semilogy(dx, error_Fouth_order_central,'c')
title('Error analysis');
xlabel('dx values');
ylabel('log_1_0 error)')
legend('Error Fouth order Right skewed','Error Fourth order Left skewed','Error Fouth order central')
function name as ''f''
function output=f(x)
output=exp(x)*cos(x);
end
post-processing analysis :
1.accuracy:
central>Left skewed >Right skewed (not always )
2.As dx values increase the accuracy is decreasing for all types.as we increasing dx values, we are ignoring higher values of `O(dx^4)` value.
let us say
`dx=.05`4th order term is about `.05^4=6.5X10^(-6)`
if use increase dx values, say dx =.1 the fourth-order term `.1^4=10^(-4)`.So we are ignoring higher value values. so error value increase.
3.its better to use log-log plot.
4. At lower values of dx, the errors are a bit high. it is due to Roundoff nature of values.it is depend on compilers and software.
let us answer this question.
In your own words, describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?
Even through accuracy wise central difference is good(Not always), but we have to encounter some problems with this method. Let see what difficulties those are.
Difficulties near the boundary:
every problem must start with some boundary conditions.
For better understanding consider one-dimensional heat conduction example.
`(d^2T)/(dx^2)=constant` say `T_1=10^oc and T_2=100^0c`
if we apply central scheme at 2nd node(fouth order approximation)
`(d^2T)/dx^2=(-30T_i+16T_(i+1)-T_(i+2)+16T_(i-1)-T_(i-2))/(12h^2)` at `i=2`
`(d^2T)/dx^2=(-30T_2+16T_(3)-T_(4)+16T_(1)-T_(0))/(12h^2)`
here `T_2,T_3,T_4` are unknow values.`T_1`is known value
What about `T_0`. is it exist?
No, physically point 'i=0' is not exist. This is the difficultly associated with the central scheme. Here mathematics is right but physically it is not possible.
so at boundaries, we have to use one-sided schemes(here right-hand skewed should use, but not left hand skewed).
Here a central scheme gives accuracy and one-sided scheme give realistic nature. Can't we use both for optimum accuracy? while we are programming.
yes, we can use. but the code will be awkward.
Thanks😊😊😊😊😊😊😊😊😊😊