Week 3 - Taylor table method and Matlab code
Objective:
- To derive the 4th order approximations of the second-order derivative by using Central difference, Skewed right-sided difference & Skewed left-sided difference and prove that your skewed schemes are fourth-order accurate.
- To write a program in Matlab evaluating the second-order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations.
Derivation of 4th order approximations of the second-order derivative
- Central difference scheme
To derive using the central difference scheme, A stencil need with a particular number of nodes needs to be taken, so that the results are of 4th order approximation.
For Central Difference Schemes
Number of nodes in stencil= p + q - 1
Where p = order of derivative
q = order of approximation
Here
P = 2
q = 4
Number of nodes = p+q-1 = 2+4-1
= 5
Therefore, the stencil of points (i-2), (i-1), i, (i+1), (i+2) is taken for the derivation of fourth-order approximation schemes.
Rewriting in Taylor table format
Converting the Taylor table in equation format and simplifying the fraction, we get
Now we have 5 equations & 5 unknowns. On converting these equations into Matrix, Matlab can solve the equations.
Rewriting this into matrix format
This matrix is of form A X = B. X can be calculated by using X = B * (Inverse of A)
Matlab code is written for solving these equations
- Skewed Right-Sided difference scheme
To derive using the Right-sided difference scheme, A stencil need with a particular number of nodes needs to be taken, so that the results are of 4th order approximation.
For Right-sided Difference Schemes
Number of nodes = p + q
Where p = order of derivative
q = order of approximation
Here
P = 2
q = 4
Number of nodes = p+q-1 = 2+4
= 6
Therfore, the sencil of points i, (i+1), (i+2), (i+3), (i+4), (i+5) is taken for derivation
Rewriting in Taylor table format
Converting the Taylor table in equation format and simplifying the fraction, we get
Now we have 5 equations & 5 unknowns. On converting these equations into Matrix, Matlab can solve the equations.
Rewriting this into matrix format
This matrix is of form A X = B. X can be calculated by using X = B * (Inverse of A)
Matlab code is written for solving these equations
- Skewed Left-Sided difference scheme
To derive using the Left-sided difference scheme, A stencil need with a particular number of nodes needs to be taken, so that the results are of 4th order approximation.
For Right-sided Difference Schemes
Number of nodes = p + q
Where p = order of derivative
q = order of approximation
Here
P = 2
q = 4
Number of nodes = p+q-1 = 2+4
= 6
Therfore, the sencil of points (i-5), (i-4), (i-3), (i-2), (i-1), i is taken for derivation
Rewriting in Taylor table format
Converting the Taylor table in equation format and simplifying the fraction, we get
Now we have 5 equations & 5 unknowns. On converting these equations into Matrix, Matlab can solve the equations.
Rewriting this into matrix format
This matrix is of form A X = B. X can be calculated by using X = B * (Inverse of A)
Matlab code is written for solving these equations
Explanations of Code - 1
function out = analytical_solution(x)
% z is considered as variable
syms z
% Equation which is to be differentiated
f = exp(z)*cos(z);
% f2 - second order differentiation of f
f2 = diff(f,2);
% substituting value of z which is to obtained from main program
out = vpa(subs(f2,z,x));
endThis code differentiates the given function to its second-order & provides the output. This output is the analytical solution for the given function.
Explanations of Code - 2
function y = f(x);
% Ddefinfng the given function
y = exp(x)*cos(x);
out = y;
end
This code evaluates for the given value of x with the function provided inside it and gives the output.
Explanations of Code - 3
clear all
close all
clc
x = pi/3;
dx = linspace(pi/4,pi/400,100);
% Central Difference
% Matrix multiplication to find values of a,b,c,d,e
% p - 5 x 5 matrix - valume from taylor table
p_c = [1 1 1 1 1; -2 -1 0 1 2; 4/2 1/2 0 1/2 4/2; -8/6 -1/6 0 1/6 8/6; 16/24 1/24 0 1/24 16/24];
% q - matrix after equal to sign
q_c = [ 0; 0; 1; 0; 0];
% r - unknown matrix
r_c = inv(p_c) * q_c ;
% Assigning values to a, b, c, d, e
a_c = r_c(1,1);
b_c = r_c(2,1);
c_c = r_c(3,1);
d_c = r_c(4,1);
e_c = r_c(5,1);
% Skewed Right - Sided Difference
% Matrix multiplication to find values of a,b,c,d,e,g
% p - 6 x 6 matrix - value from taylor table
p_r = [1 1 1 1 1 1; 0 1 2 3 4 5; 0 1/2 4/2 9/2 16/2 25/2; 0 1/6 8/6 27/6 64/6 125/6; 0 1/24 16/24 81/24 256/24 625/24; 0 1/120 32/120 243/120 1024/120 3125/120];
% q - matrix after equal to sign
q_r = [ 0; 0; 1; 0; 0; 0];
% r - unknown matrix
r_r = inv(p_r)* q_r;
% Assigning values to a, b, c, d, e, g
a_r = r_r(1,1);
b_r = r_r(2,1);
c_r = r_r(3,1);
d_r = r_r(4,1);
e_r = r_r(5,1);
g_r = r_r(6,1);
% Skewed Left - Sided Difference
% Matrix multiplication to find values of a,b,c,d,e,g
% p - 6 x 6 matrix - value from taylor table
p_l = [1 1 1 1 1 1; -5 -4 -3 -2 -1 0; 25/2 16/2 9/2 4/2 1/2 0; -125/6 -64/6 -27/6 -8/6 -1/6 0; 625/24 256/24 81/24 16/24 1/24 0; -3125/120 -1024/120 -243/120 -32/120 -1/120 0];
% q - matrix after equal to sign
q_l = [ 0; 0; 1; 0; 0; 0];
% r - unknown matrix
r_l = inv(p_l)* q_l;
% Assigning values to a, b, c, d, e
a_l = r_l(1,1);
b_l = r_l(2,1);
c_l = r_l(3,1);
d_l = r_l(4,1);
e_l = r_l(5,1);
g_l = r_l(6,1);
% preallocating for faster output
CDS_solution = zeros(1,100);
CDS_error = zeros(1,100);
RSS_solution = zeros(1,100);
RSS_error = zeros(1,100);
LSS_solution = zeros(1,100);
LSS_error = zeros(1,100);
a = zeros(1,100);
b = zeros(1,100);
c = zeros(1,100);
d = zeros(1,100);
e = zeros(1,100);
h = zeros(1,100);
j = zeros(1,100);
k = zeros(1,100);
l = zeros(1,100);
m = zeros(1,100);
for i=1 : length(dx)
% Allocating Taylor series nodes at various points
a(i)= f(x - (5 * dx(i)));
b(i)= f(x - (4 * dx(i)));
c(i)= f(x - (3 * dx(i)));
d(i)= f(x - (2 * dx(i)));
e(i)= f(x - (1 * dx(i)));
h(i)= f(x + (1 * dx(i)));
j(i)= f(x + (2 * dx(i)));
k(i)= f(x + (3 * dx(i)));
l(i)= f(x + (4 * dx(i)));
m(i)= f(x + (5 * dx(i)));
% CDS Error
CDS_solution(i) = ((a_c * d(i)) + (b_c * e(i)) + (c_c * f(x)) + (d_c *h(i)) + (e_c * j(i)))/ (dx(i)^2);
CDS_error(i) = abs(CDS_solution(i) - analytical_solution (x));
% Right side skewed approximation Error
RSS_solution(i) = ((a_r * f(x)) + (b_r * h(i)) + (c_r * j(i)) + (d_r * k(i)) + (e_r * l(i)) + (g_r * m(i)))/ (dx(i)^2);
RSS_error(i) = abs(RSS_solution(i) - analytical_solution (x));
% Left side skewed approximation Error
LSS_solution(i) = ((a_l * a(i)) + (b_l * b(i)) + (c_l * c(i)) + (d_l * d(i)) + (e_l * e(i)) + (g_l * f(x)))/ (dx(i)^2);
LSS_error(i) = abs(LSS_solution(i) - analytical_solution (x));
end
figure(1)
plot(dx,CDS_error,'r')
hold on
plot(dx,RSS_error,'g')
hold on
plot(dx,LSS_error,'b')
grid on
xlabel('dx')
ylabel('Absolute error')
title ('Consistency of Different Schemes')
legend('central difference approximation error','skewed right sided approximation error','skewed left sided approximation error')
The main program starts with default commands on clearing all the stored values & the screen and closing the windows.
The value of x & dx is assigned initially.
Then the values on the Central Difference Scheme are calculated.
- Values from Taylor Table are given in matrix format followed by the summation value in the table.
- The unknown values are calculated inside the Matlab by multiplying the inverse matrix of values from the Taylor table to the summation value.
- Now each value of the unknown matrix is assigned to separate variables that are to be used in the program.
Similarly, the unknown values of the Right side sked difference scheme & Left Side Skewwd difference Schemes are calculated and assigned to separate variables.
As a suggestion from Matlab, for increasing the computing speed, the values such as CDS_solution, CDS_Error & similar variables are preallocated with the size of the matrix.
Inside a for loop, for every value of dx
- The value of the function (given in code 2) is called and its values are assigned to separate variables.
- The value of the Central Difference Scheme solution, Central Difference Scheme error, Right side skewed difference scheme solution, Right side skewed difference scheme error, Left Side Skewed difference Schemes solution, Left Side Skewed difference Schemes error is calculated.
All the values of Central Difference Scheme error, Right side sked difference scheme error, Left Side Skewed difference Schemes error are plotted against dx inside the same figure. followed by suitable Labels, Titles & Legend for the plot.
Results & Discussion
From the plot, we infer that, within the specified limit, the central difference scheme is more accurate when compared with the skewed schemes. When comparing with the left & right side skewed system, Left side skewed stencils produce more accurate results whereas the produce accurate result for a limited value range of dx and following it, the error values increases exponentially.
What a skewed scheme can do
In a practical aspect, For computing the values at a point in the middle of the model or when computing from the boundary conditions, It is impossible to go with central difference schemes as the values of after nodes are not known. Hence in that condition, only skewed stencils can be used to find the values of adjacent nodes. This is the benefit skewed schemes over central difference schemes.