Discretization basics for range of dx

OBJECTIVE: To write a program that compares the first, second and fourth order approximations of the first derivative against the analytical or exact derivative with a varying value of dx.

The function that is to be computed is f(x) = sin(x)/x^3 and the given value for x is x=pi/3.

Numerical computation for first order approximation:

forward differencing = `(f(x+a)-f(x))/dx`

Numerical computation for second order approximation:

central differencing = `(f(x+a)-f(x-a))/(2*dx)`

Numerical computation for fourth order approximation:

central differencing = `(f(x-2*dx)-(8*f(x-dx))+(8*f(x+dx))-f(x-2*dx))/(12*dx)`

Now we define a function for first order approximation. Here we calculate the analytical derivative value and the first order approximation value using forward differencing. Then an error term is introduced which is the absolute difference between forward differencing value and the analytical value.

FIRST ORDER APPROXIMATION FUNCTION CODE:

% Defining the function for the first order approximation

function error_first_order = first_order_approximation(x,dx)
  
% analytical function = sin(x)/x^3;

% analytical derivative

% f'(x) = (x^3*(cos(x))-sin(x)*3*x^2)/x^6;


analytical_derivative = (x^3*(cos(x))-sin(x)*3*x^2)/x^6;


%first order approximation

forward_differencing = ((sin(x+dx)/(x+dx)^3) - (sin(x)/(x)^3))/dx;

% Storing the output values

error_first_order = abs(forward_differencing - analytical_derivative)

end

Similarly we define a function for second order approximation and here the error is the absolute difference between central differencing value and the analytical derivative value.

SECOND ORDER APPROXIMATION FUNCTION CODE:

% Defining the function for the second order approximation

function error_second_order = second_order_approximation(x,dx)
  
% analytical function = sin(x)/x^3;

% analytical derivative

% f'(x) = (x^3*(cos(x))-sin(x)*3*x^2)/x^6;

analytical_derivative = (x^3*(cos(x))-sin(x)*3*x^2)/x^6;

%second order approximation

central_differencing = ((sin(x+dx)/(x+dx)^3) - (sin(x-dx)/(x-dx)^3))/(2*dx);

% Storing the output values

error_second_order = abs(central_differencing - analytical_derivative)

end

Then we proceed to define a function for fourth order approximation and here also central differencing method is used as like second order approximation. The absolute difference between central differencing value and analytical derivative value gives us the error.

FOURTH ORDER FUNCTION CODE:

% Defining the function for the fourth order approximation

function out_fourth_order = fourth_order_approximation(x,dx)
  
% analytical function = sin(x)/x^3;

% analytical derivative

% f'(x) = (x^3*(cos(x))-sin(x)*3*x^2)/x^6;

analytical_derivative = (x^3*(cos(x))-sin(x)*3*x^2)/x^6;

%fourth order approximation

central_differencing_1 = ((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);

% Storing the output values

out_fourth_order = abs(central_differencing_1 - analytical_derivative)

end

After defining all the functions we now proceed for the main program.

At x=pi/3,we have to calculate the analytical derivative of the given function by differentiating once.

Now we have to declare the term 'dx'. Since it is been asked to compute the approximations for varying 'dx' value, we use the linspace command and declare 20 values ranging from pi/4 to pi/400.

The approximation values keep on varying as 'dx' changes and hence a for loop is introduced for the iteration process. In this loop, an iterative variable 'i' is declared and the iterative process is as follows:-

1. when i=1 then dx value is dx(1) = pi/4. so for x = pi/3 and dx = pi/4 the first order approximation value is computed and stored in first order error (1).
2. Now i gets incremented to 2 i.e. i=2 and dx becomes dx(2). With this first order error (2) value is calculated and the process continues.
3. for i=20 , dx(20) = pi/400 and first order error (20) value is obtained.
4. similarly second and fourth order errors are also calculated.

After the error values are obtained, they have to be plotted using the loglog plot for the sake of comparison regarding which order of approximation gives the least error.

MAIN PROGRAM:

  clear all
  close all
  clc
  
  % Declaring the x value
  
  x = pi/3;
  
  % Declaring the dx value with a range 
  
  dx = linspace(pi/4,pi/400,20);
  
  % To begin a for loop for iteration process and declaring variable 'i'
  
  for i  = 1:length(dx)
    
    first_order_error(i) = first_order_approximation(x,dx(i));   % calling the first order approxiamtion function
    
    second_order_error(i) = second_order_approximation(x,dx(i));  % calling the second order approxiamtion function
    
    fourth_order_error(i) = fourth_order_approximation(x,dx(i));  % calling the fourth order approxiamtion function
    
  end
  
  % Plotting the above obtained results
  
  figure(1)
  
  % using loglog command and specifying colors for the plot
  
  loglog(dx,first_order_error,'b','linewidth',3)
  
  hold on
  
  loglog(dx,second_order_error,'r','linewidth',3)
  
  hold on
  
  loglog(dx,fourth_order_error,'g','linewidth',3)
  
  legend('first order error','second order error','fourth order error')

LOG-LOG PLOT REPRESENTATION:

CONCLUSIONS:

1. From the above plot representation we can conclude that fourth order error is the least when compared to first and second order errors.
2. Though for initial values of 'dx' there is an increase in error in fourth order but  as 'dx' value decreases error also a decreases.

slope calculations for the above plot: Let us calculate slope for the initial and final points of the plot.

  i) For first order error: Point 1 = (x1,y1) = (pi/4,0.96466)   where y1 = first order error at                                                                                                                                      dx =  pi/4.     

                                   Point 2 = (x2,y2) = (pi/400,0.019454)       where y2 = first order error at                                                                                                                                   dx = pi/400.

                                   Slope = (y2-y1)/(x2-x1) = 1.2156

ii)  For second order:  Point 1 = (x1,y1) = (pi/4,7.3578)     where y1 = second order error at                                                                                                                                  dx =  pi/4.

                                 Point 2 = (x2,y2) = (pi/400,0.00019599)    where y2 = second order error at                                                                                                                           dx = pi/400.

                                 Slope = (y2-y1)/(x2-x1) = 9.4628

iii)  For fourth order:    Point 1 = (x1,y1) = (pi/4,10.019)         where y1 = fourth order error at                                                                                                                                  dx = pi/4.              

                                  Point 2 = (x2,y2) = (pi/400,0.000000066150)  where y2 = fourth order error                                                                                                                           at dx = pi/400.

                                    Slope = (y2-y1)/(x2-x1) = 12.885

ERROR VALUES:

                  Error values at dx = pi/4. 

               Error values at dx = pi/400.

From the above slope calculations we can say that slope of fourth order error is higher than second and first order error. Since slope is high therefore error is least.