Numerical discretization of various differencing schemes and error calculation for a range of dx

Objectives:

To compute the first,second and fourth order approximations for a range of dx values and plot the error vs dx.

Given:

To compute the numerical approximations of the function f(x) =`sinx/x^3` at x = pi/3 for the range of dx values.

The numerical derivatives can be computed with the help of taylor's series expansion formula

The first,second and fourth order approximations can be calculated using one of the methods

• Forward differencing method
• Backward differencing method
• Central differencing method

First order approximation - Forward differencing method

Second order approximation- Central differencing method

Fourth order approximation - Central differencing method

Formula:

Forward differencing >> `f'(x)=[f(x+dx)-f(x)]/dx`

Backward differencing >> `f'(x)=[f(x)-f(x-dx)]/dx`

Central differencing >>`f'(x)=[f(x+dx)-f(x-dx)]/(2dx)`

Octave Program:

Function to calculate first order approximation

%Function to calculate the first_order_approximation
function result=first_order_approximation(x,dx)
  
  exact_derivative=((cos(x)*x^3)-(sin(x)*3*(x^2)))/x^6;
%forward_differencing
%first_order_approximation f'(x)=(f(x+dx)-f(x))/dx

  first_order=((sin(x+dx)/(x+dx)^3)-(sin(x)/x^3))/dx;
  result=abs(first_order-exact_derivative);
end

Function to calculate second order approximation

%Function to calculate the second_order_approximation
function result=second_order_approximation(x,dx)
  exact_derivative=((cos(x)*x^3)-(sin(x)*3*(x^2)))/x^6;
%Central_differencing
%second_order_approximation f'(x)=(f(x+dx)-f(x-dx))/2dx
  second_order=((sin(x+dx)/(x+dx)^3)-(sin(x-dx)/(x-dx)^3))/(2*dx);
  result=abs(second_order-exact_derivative);
end

Function to calculate fourth order approximation

%Function to calculate the fourth_order_approximation
function result=fourth_order_approximation(x,dx)
  
  exact_derivative=((cos(x)*x^3)-(sin(x)*3*(x^2)))/x^6;
%Central_differencing
%first_order_approximation f'(x)=(f(x+dx)-f(x))/dx

  fourth_order=((sin(x-2*dx)/(x-2*dx)^3)-(8*(sin(x-dx)/(x-dx)^3))+(8*(sin(x+dx)/(x+dx)^3))-(sin(x+2*dx)/(x+2*dx)^3))/(12*dx);
  result=abs(fourth_order-exact_derivative);
endfunction

Main code:

clear all;
close all;
clc;
%The function f(x)=sin(x)/x^3 is computed at x=pi/3
x=pi/3;
%dx is an array of values to compute the error of numerical scheme for a range of dx values
dx=linspace(pi/4,pi/4000,20);
%Looping through the range of dx values to calculate the error
for i =1:length(dx)
  first_order_error(i)=first_order_approximation(x,dx(i));
  second_order_error(i)=second_order_approximation(x,dx(i));
  fourth_order_error(i)=fourth_order_approximation(x,dx(i));
end

%Plotting the values for the error vs dx to compare the numerical schemes
figure(1)
plot(dx,first_order_error,'r')
hold on
plot(dx,second_order_error,'b')
hold on
plot(dx,fourth_order_error,'g')
legend('first order error','second order error','fourth order error');
legend('location','southeast');
xlabel('dx');
ylabel('error')

%Log plot is used as there is difference in the error terms of numerical schemes
figure(2)
loglog(dx,first_order_error,'r')
hold on
loglog(dx,second_order_error,'b')
hold on
loglog(dx,fourth_order_error,'g')
legend('first order error','second order error','fourth order error');
legend('location','southeast');
xlabel('dx');
ylabel('error')

%Slope is calculated for the curves at all points
slope1=diff(first_order_error)./diff(dx);
slope2=diff(second_order_error)./diff(dx);
slope4=diff(fourth_order_error)./diff(dx);

Variation of the error for the range of dx:

From the above plot we cannot identify the variation for the error for the various numerical schemes since first order error is in range of `e^-1` where as fourth order error is range of `e^-12` so there is no variation and hence log plot is used.

Log plot :

From the above log plot we can clearly see that error of first order approximation is higher than the second and fourth order apppriximations.The error of fourther order approximation reaches to zero as the dx value is reduced(grid space is reduced).

In fourth order approximation the absolute error is low as dx value gets reduced, the error is maximum at higher dx values.Fourth order approximation reaches zero at rapid rate than other order approximations. In order to get the approximate solution , higher order approximation is preffered.

Fourth order approximation>Second order approximation>First order approximation.

Conclusion:

• Fourth order approximation reaches zero error at a rapid rate than the other numerical approximations for range of dx.
• dx is an important term in numerical discretization to attain the approximation solution.
• dx must be selected correctly to get the approximate solution, or else it results in divergence of the error.
• Higher order approximation is preffered as the approximate solution is obtained for the given function, it may vary based on the function.
• In 3-D if the mesh is finer (dx is small ) the accuracy of the result is high , or else the result diverges (dx is large).