Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
CENTRE DIFFERENCE SCHEME-
So, for central difference scheme formula for choosing no nodes is
No of nodes = order of approximation + degree of derivative -1
here, order of approximation = 4 ; order of derivative = 2
hence, no of nodes = 4 + 2 -1
= 5
Now in a central difference scheme there should be equal no of nodes on left and right hand side as shown in above figure
`(del^2f)/(delx^2)`= af(i-2) + bf(i-1) + cf(i) +df(i+1) + ef(i+2) [ a,b,c,d,e are unknown coefficients]
Now expanding each of the derivative using taylor's series:
a*f(i-2) = a[f(x) - (`2trianglex`f'(x))/1! + `(2trianglex)^2`f''(x)/2! - `(2trianglex)^3`f'''(x)/3! + `(2trianglex)^4`f''''(x)/4! - `(2trianglex)^5`f'''''(x)/5! + `(2trianglex)^6`f''''''(x)/6!]
b*f(i-1) = b[f(x) - (`trianglex`f'(x))/1! + `(trianglex)^2`f''(x)/2! - `(trianglex)^3`f'''(x)/3! + `(trianglex)^4`f''''(x)/4! - `(trianglex)^5`f'''''(x)/5! + `(trianglex)^6`f''''''(x)/6!]
c*f(i) = c[f(i)]
d*f(i+1) = d[f(x) + (`trianglex`f'(x))/1! + `(trianglex)^2`f''(x)/2! + `(trianglex)^3`f'''(x)/3! + `(trianglex)^4`f''''(x)/4! + `(trianglex)^5`f'''''(x)/5! + `(trianglex)^6`f''''''(x)/6!]
e*f(i+2) = e[f(x) + (`2trianglex`f'(x))/1! + `(2trianglex)^2`f''(x)/2! + `(2trianglex)^3`f'''(x)/3! + `(2trianglex)^4`f''''(x)/4! + `(2trianglex)^5`f'''''(x)/5! + `(2trianglex)^6`f''''''(x)/6!]
we can also assume x=i , just to stop any cofusion between x & `trianglex`
now addning up the coefficient
a + b+ c+ d+ e = 0 ..............1
-2a -b +d +2e = 0 ..................2
4a/2 + b/2 + d/2 + 4e/2 = 1 .............3
-8a/6 - b/6 + d/6 + 8e/6 = 0 ..................4
16a/24 + b/24 + d/24 + 16e/24 = 0 ...........5
use the above equation for finding out the coefficient a,b,c,d,e
form a matrix of above equation first
A=
B = [a;b;c;d;e]
By performing inverse of matrix AB = X in matlab we can find out value for each coefficient
use this inbuilt matlab function B = mldivide(A,X)
a = -0.0833, b = 1.3333 , c= -2.5 , d = 1.3333, e = -0.0833
here a = e & b = d
putting value of coefficient in f'''''(i) & f''''''(i) to find out value of ? & ?
f'''''(i) = (-32a/120 -b/120 + d/120 32e/120)
= 0 ( sum of its coefficient comes to 0)
f''''''(i) = (64a/720 + b/720 + d/720 + 64e/720) here, after putting values of coefficient our equation doesn't equals to zero, Hence sum of its coefficient = 1
Only for f''(i) & f''''''(i), sum of coefficient = 1 for rest it is 0
Now adding f''(i) & f''''''(i)
`(del^2f)/(delx^2)`= `(trianglex^2)`f''(i) + `(trianglex)^6`f''''''(i)
= (2a + b/2 + d/2 + 2e)`(trianglex^2)`f''(i) + (64a/720 + b/720 + d/720 + 64e/720)`(trianglex)^6`f''''''(i)
since, `(del^2f)/(delx^2)` = af(i-2) + bf(i-1) + cf(i) +df(i+1) + ef(i+2)
hence, af(i-2) + bf(i-1) + cf(i) +df(i+1) + ef(i+2) = (2a + b/2 + d/2 + 2e)`(trianglex)^2`f''(i) + (64a/720 + b/720 + d/720 + 64e/720)`(trianglex)^6`f''''''(i)
Dividing LHS and RHS by `trianglex^2`
we get,
(af(i-2) + bf(i-1) + cf(i) +df(i+1) + ef(i+2))/`(trianglex)^2`= (2a + b/2 + d/2 + 2e)f''(i) + (64a/720 + b/720 + d/720 + 64e/720)`(trianglex)^4`f''''''(i)
(af(i-2) + bf(i-1) + cf(i) +df(i+1) + ef(i+2))/`(trianglex)^2` = 0.9997f''(i) + (-0.0111)`(trianglex)^4`f''''''(i)
LEFT HAND SKEWED SCHEME-
So, for right hand scheme formula for choosing no nodes is
No of nodes = order of approximation + degree of derivative -1
here, order of approximation = 4 ; order of derivative = 2
hence, no of nodes = 4 + 2
= 6
Now in a central difference scheme there should be equal no of nodes on left and right hand side as shown in above figure
`(del^2f)/(delx^2)` = af(i-5) + bf(i-4) + cf(i-3) + df(i-2) + ef(i-1) + gf(i)
af(i-5) = a[f(i) - `(5trianglex)`f'(x)/1! + `(5trianglex)^2`f''(i)/2! - `(5trianglex)^3`f'''(i)/3! + `(5trianglex)^4`f''''(i)/4! - `(5trianglex)^5`f'''''(i)/5! + `(5trianglex)^6`f''''''/6!]
bf(i-4) = b[f(i) - `(4trianglex)`f'(i)/1! + `(4trianglex)^2`f''(i)/2! - `(4trianglex)^3`f'''(i)/3! -`(4trianglex)^4`f''''(i)/4! - `(4trianglex)^5`f'''''(i)/5! + `(4trianglex)^6`f''''''(i)/6!]
cf(i-3) = c[f(i) - `(3trianglex)`f'(i)/1! + `(3trianglex)^2`f''(i)/2! - `(3trianglex)^3`f'''(i)/3! + `(3trianglex)^4`f''''(i)/4! - `(3trianglex)^5`f'''''(i)/5! + `(3trianglex)^6`f''''''(i)/6!]
df(i-2) = d[f(i) - `(2trianglex)`f'(i)/1! + `(2trianglex)^2`f''(i)/2! - `(2trianglex)^3`f'''(i)/3! + `(2trianglex)^4`f''''(i)/4! - `(2trianglex)^5`f'''''(i)/5! + `(2trianglex)^6`f''''''(i)/6!]
ef(i-1) = e[f(i) - `(trianglex)`f'(i)/1! + `(trianglex)^2`f''(i)/2! - `(trianglex)^3`f'''(i)/3! + `(trianglex)^4`f''''(i)/4! - `(trianglex)^5`f''''(i)/5! + `(trianglex)^6`f''''''(i)/6!]
gf(i) = gf(i)
Now adding up the coefficients
a + b + c +d +e + g = 0 .....................1
-5a -4b -3c -2b -e = 0 ......................2
25a/2 + 16b/2 + 9c/2 + 4d/2 + 1e/2 = 1 ...........3
-125a/6 - 64b/6 - 27c/6 - 8d/6 -e/6 = 0 ............4
5^4a/24 + 4^4b/24 + 91c/24 + 16d/24 + e/24 = 0..............5
-5^5a/120 - 4^5b/120 - 3^5c/120 - 2^5d/120 - e/120 = 0.........6
By performing inverse of matrix AB = X in matlab we can find out value for each coefficient
use this inbuilt matlab function B = mldivide(A,X)
Value of coefficients are
a= -0.8333
b= 5.0833
c= -13.000
d = 17.8330
e= -12.8333
g= 3.7500
putting these values in sum of coefficient of f''''''(i) we get a non zero value
hence, f''''''(i) = 1
sum of coefficient of f''(i) & f''''''(i) is 1 and sum of rest sre zero
NOW, `(delf^2)/(del^2x)` = (25a/2 + 16b/2 + 9c/2 + 4d/2 + e/2)`(trianglex)^2`f''(i) + [(5^6)a/720+ (4^6)b/720 + (3^6)c/720 + (2^6)d/720 + e/720]`(trianglex)^6`f''''''(i)
putting values of a,b,c,d,e
`(delf^2)/(del^2x)` = (25*-0.0833/2 + 16*5.0833/2 + 9*-13.000/2 + 4*17.8330/2 + -12.8333/2)`(trianglex)^2`f''(i)+ [(5^6)*-0.8333/720+ (4^6)*5.08333/720 + (3^6)*-13.000/720 + (2^6)*12.8333/720 + -12.8333/720]`(trianglex)^6`f''''''(i)
since `(delf^2)/(del^2x)` = af(i-5) + bf(i-4) + f(i-3) + f(i-2) + f(i-1) + f(i)
hence, af(i-5) + bf(i-4) + f(i-3) + f(i-2) + f(i-1) + f(i) = (25*-0.0833/2 + 16*5.0833/2 + 9*-13.000/2 + 4*17.8330/2 + -12.8333/2)`(trianglex)^2`f''(i) + [(5^6)*-0.8333/720+ (4^6)*5.08333/720 + (3^6)*-13.000/720 + (2^6)*12.8333/720 + -12.8333/720]`(trianglex)^6`f''''''(i)
divinding by `(trianglex)^2`
[af(i-5) + bf(i-4) + f(i-3) + f(i-2) + f(i-1) + f(i)]/`(trianglex)^2` = (25*-0.0833/2 + 16*5.0833/2 + 9*-13.000/2 + 4*17.8330/2 + -12.8333/2)f''(i) + [(5^6)*-0.8333/720+ (4^6)*5.08333/720 + (3^6)*-13.000/720 + (2^6)*12.8333/720 + -12.8333/720]`(trianglex)^4`f''''''(i)
PDE = AE , partial differential equation is now converted to algebric equation
RIGHT HAND SKEWED SCHEME-
So, for right hand scheme formula for choosing no nodes is
No of nodes = order of approximation + degree of derivative -1
here, order of approximation = 4 ; order of derivative = 2
hence, no of nodes = 4 + 2
= 6
Now in a central difference scheme there should be equal no of nodes on left and right hand side as shown in above figure
`(delf^2)/(delx^2)` = af(i) + bf(i+1) + cf(i+2) +df(i+3)+ef(i+4) + gf(i+5)
af(i) = af(i)
bf(i+1) = b[f(i) + `(trianglex)`f'(i)/1! + `(trianglex)^2`f''(i)/2! + `(trianglex)^3`f'''(i)/3! + `(trianglex)^4`f''''(i)/4! + `(trianglex)^5`f''''(i)/5! + `(trianglex)^6`f''''''(i)/6!]
cf(i+2) = c[f(i) + `(2trianglex)`f'(i)/1! + `(2trianglex)^2`f''(i)/2! + `(2trianglex)^3`f'''(i)/3! + `(2trianglex)^4`f''''(i)/4! + `(2trianglex)^5`f'''''(i)/5! + `(2trianglex)^6`f''''''(i)/6!]
df(i+3) = d[f(i) + `(3trianglex)`f'(i)/1! + `(3trianglex)^2`f''(i)/2! + `(3trianglex)^3`f'''(i)/3! + `(3trianglex)^4`f''''(i)/4! + `(3trianglex)^5`f'''''(i)/5! + `(3trianglex)^6`f''''''(i)/6!]
ef(i+4) = e[f(i) + `(4trianglex)`f'(i)/1! + `(4trianglex)^2`f''(i)/2! + `(4trianglex)^3`f'''(i)/3! +`(4trianglex)^4`f''''(i)/4! + `(4trianglex)^5`f'''''(i)/5! + `(4trianglex)^6`f''''''(i)/6!]
gf(i+5) = g[f(i) + `(5trianglex)`f'(x)/1! + `(5trianglex)^2`f''(i)/2! + `(5trianglex)^3`f'''(i)/3! + `(5trianglex)^4`f''''(i)/4! + `(5trianglex)^5`f'''''(i)/5! + `(5trianglex)^6`f''''''/6!]
now adding up the coefficient
a + b +c +d+e+g = 0..............1
b + 2c +3d +4e +5g = 0 ............2
b/2 + 4c/2 + 9d/2 + 16e/2 + 25g/2 = 1 .....3
b/6 + 8c/6 + 27d/2 + 64e/6 + 125g/6 = 0 ...............4
b/24 + 16c/24 + 91d/24 + (4^4)e/24 + (5^4)g/24= 0 ..............5
b/120 + 32c/120 + (3^5)d/120 + (4^5)e/120 + (5^5)g/120 = 0..................6
By performing inverse of matrix AB = X in matlab we can find out value for each coefficient
use this inbuilt matlab function B = mldivide(A,X)
hence the value of a,b,c,d,e,g are
a=3.7500
b=-12.8333
c=17.8333
d=-13.000
e=5.08333
g=-0.8333
putting these values in sum of coefficient of f''''''(i) we get a non zero value
hence, f''''''(i) = 1
sum of coefficient of f''(i) & f''''''(i) is 1 and sum of rest are zero
NOW, `(delf^2)/(del^2x)` = (b/2 +4c/2 + 9d/2 +16e/2 + 25g/2) `(trianglex)^2`f''(i) + [b/720+ 64c/720 + (3^6)d/720 + (4^6)e/720 + (5^6)g/720]`(trianglex)^6`f''''''(i)
putting values of a,b,c,d,e
`(delf^2)/(del^2x)` = (-12.8333/2 +4*17.8333/2 + 9*-13.000/2 +16*5.08333/2 + 25*-0.8333/2) `(trianglex)^2`f''(i) + [-12.8333/720+ 64*17.8333/720 + (3^6)*-13.000/720 + (4^6)*5.08333/720 + (5^6)*-0.8333/720 `(trianglex)^6`f''''''(i)
since `(delf^2)/(del^2x)` = af(i-5) + bf(i-4) + f(i-3) + f(i-2) + f(i-1) + f(i)
hence, af(i-5) + bf(i-4) + f(i-3) + f(i-2) + f(i-1) + f(i) = `(trianglex)^2`f''(i) +
`(trianglex)^6`f''''''(i)
divinding by `(trianglex)^2`
[af(i-5) + bf(i-4) + f(i-3) + f(i-2) + f(i-1) + f(i)]/`(trianglex)^2` = (-12.8333/2 +4*17.8333/2 + 9*-13.000/2 +16*5.08333/2 + 25*-0.8333/2)f''(i) + (-12.8333/720+ 64*17.8333/720 + (3^6)*-13.000/720 + (4^6)*5.08333/720 + (5^6)*-0.8333/720 )`(trianglex)^4`f''''''(i)
PDE = AE , partial differential equation is now converted to algebric equation
MATLAB CODE FOR FOURTH ORDER APPROXIMATION OF 2ND ORDER DERIVATIVE-
close all
clear all
clc
%analytical function
%f(x) = exp(x)*cos(x)
%derivative of this function
%f''(x) = -2exp(x)*sin(x)
x= pi/3;
dx = linspace(pi/3,pi/3000,30);
analytical_derivative = -2*exp(x)*sin(x);
%for central differencing value of coefficients
a = -0.0833333333333;
b = 1.333333333333;
c = -2.5;
d = 1.333333333333;
e = -0.0833333333333;
%for left hand skewed scheme
a1 = -0.83333333333;
b1 = 5.083333333333;
c1 = -13.0000;
d1 = 17.83333333333;
e1 = -12.83333333333;
g1 = 3.75000;
%for right hand skewed scheme
a2 = 3.75000;
b2 = -12.83333333333;
c2 = 17.83333333333;
d2 = -13.0000;
e2 = 5.083333333333;
g2 = -0.83333333333;
for i = 1:length(dx)
Dx = dx(i);
central_difference_method(i) = (a.*(exp(x-2*Dx).*cos(x-2*Dx)) + b.*(exp(x-1*Dx).*cos(x-1*Dx)) + c*(exp(x)*cos(x))+ d.*(exp(x+1*Dx).*cos(x+1*Dx))+ e.*(exp(x+2*Dx).*cos(x+2*Dx)))./(Dx.^2);
left_hand_skewed_scheme(i) = (a1.*(exp(x-5*Dx).*cos(x-5*Dx)) + b1.*(exp(x-4*Dx).*cos(x-4*Dx)) + c1.*(exp(x-3*Dx).*cos(x-3*Dx)) + d1.*(exp(x-2*Dx).*cos(x-2*Dx)) + e1.*(exp(x-1*Dx).*cos(x-1*Dx)) + g1.*(exp(x).*cos(x)))./(Dx.^2);
right_hand_skewed_scheme(i) = (a2*(exp(x)*cos(x)) + b2.*(exp(x+1*Dx).*cos(x+1*Dx)) + c2.*(exp(x+2*Dx).*cos(x+2*Dx)) +d2.*(exp(x+3*Dx).*cos(x+3*Dx))+e2.*(exp(x+4*Dx).*cos(x+4*Dx)) + g2.*(exp(x+5*Dx).*cos(x+5*Dx)))./(Dx.^2);
end
%calculation for errors
central_difference_error = abs(analytical_derivative - central_difference_method);
left_hand_error = abs(analytical_derivative - left_hand_skewed_scheme);
right_hand_error = abs(analytical_derivative - right_hand_skewed_scheme);
%plotting of dx Vs errors
loglog(dx,central_difference_error,'linewidth',3,'color','r')
hold on
loglog(dx,left_hand_error,'linewidth',3,'color','b')
hold on
loglog(dx,right_hand_error,'linewidth',3,'color','g')
grid on
xlabel('dx')
ylabel('truncation error')
legend('central diff error','left hand error','right hand error','location','southeast')
Before starting to explain this i would to inform that there are 2 methods by which you can solve this code
(i) by using and creating function for each scheme
(ii) by using loop
We are using here (ii) method as it is less lengthy and much easier to write and understand, In loop method we will run a loop in which we will cover all there schemes simultaneously and we dont have write seperate functions for each scheme and then composing each function into a single to find out the output errors.
Explaination of code:
- First specify the function which is exp(x)cos(x)
- Find out 2nd order derivative of this function, which means differentiate it two times. we get f''(x) = -2exp(x)sin(x)
- specify the value for x and dx
- use linspace command to choose 30 values of dx between pi/30 and pi/300
- Now first find value analytical derivative of function by putting value of x in 2nd order derivative
- specify the values of coefficient for all three schemes
- Now we will apply loop and inside of this loop we will write equation for central, righ hand & left hand scheme
- We will use the equation which we have derived in our schemes which is explained in detail above, kindly go through it , it will help you in under standing the code
- Algebric equation for central differnce scheme = af(i-2) + bf(i-1) + cf(i) +df(i+1) + ef(i+2)
- Algeric equation for left hand skewed scheme = af(i-5) + bf(i-4) + cf(i-3) + df(i-2) + ef(i-1) + gf(i)
- Algebric equation for right hand skewed scheme = af(i) + bf(i+1) + cf(i+2) +df(i+3)+ef(i+4) + gf(i+5)
- f is function f(x) = exp(x)cos(x), f(i-2) is written as exp(x-2*dx)*cos(x-2*dx), (i = x)
- for values of coefficients increase no of places upto 10 points after decimals, incrase value of coefficient after decimals in order to get better result, ihave used 10 points after decimal
- loop is used as it will execute the value of equation at each value of dx by taking one value of dx at a time
- put Dx = dx(i), so that we dont have to use dx(i) at each place in our equations
- write (i) after central_difference_scheme so that all 30 values of central_difference scheme are displayed in workspace window
- calculate the values for errors for all three schems.
- use abs command , which is an inbuilt command in matlab to find out absolute difference value between analytical approach and approximation approach
- now we have to plot the values of errors we got Vs dx
- use loglog command to plot dx Vs errors
- use hold on command, this command is used when we have to ploy more than one plot in a single figure, this command avoids over writing of new plot over previous plot and help to keep all plot in single figure.
- provide labels and legends, legends are used to give identity to each plot in a figure
- in legend command 'location' is used to select position of legend dialogue box in the figure.
Graph for our code:
