MATLAB Program - Numerical Discretization for a Range of dx

Aim: Write a program that compares first order, second order and fourth order derivatives.

Objectives:

• To compute analytical derivative of given function.
• Compare this to numerical methods and evaluate the absolute error in various schemes

Theory:

In applied mathematics, discretization is the process of transferring continuous functions, models, variables, and equations into discrete counterparts. This process is usually carried out as a first step towards making them suitable for evaluation and/or implemetation in computing techniques.

In this program, we'll compare first, second and fourth order approximations to the analytical solution.

Given Function: `f(x) = (sin(x))/(x^3)`

Derivative at: `x = (pi)/3`

Code:

The program is as shown below. It is a MATLAB program.

1. Main Program:

%Discretization Basics
clear all
close all
clc

x = pi/3;
dx = linspace(pi/4,pi/4000,100);

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));
endfor

figure(1)
hold on
grid on
loglog(dx,first_order_error)
loglog(dx,second_order_error)
loglog(dx,fourth_order_error)
legendstr = {'first order error','second order error','fourth order error'};
legend(legendstr)
xlabel('dx')
ylabel('error')

2. Functions for calculating approximations

function out = first_order_approximation(x,dx)
  %given function f(x) = sin(x)/x^3
  analytical_derivative = (x^3*cos(x)-sin(x)*3*x^2)/x^6;

  %First order approximation using forward differencing scheme
  % forward_differincing = (f(x+dx) - f(x))/dx 
  forward_differencing = ((sin(x+dx))/((x+dx)^3) - (sin(x))/(x^3))/dx;
  out = abs(forward_differencing - analytical_derivative);
endfunction

function out = second_order_approximation(x,dx)
  %given function f(x) = sin(x)/x^3
  analytical_derivative = (x^3*cos(x)-sin(x)*3*x^2)/x^6;

  %Second order approximation using central differencing scheme.
  % central_differincing = (f(x+dx) - f(x-dx))/2*dx 
  central_differencing_second_order = ((sin(x+dx))/((x+dx)^3)-(sin(x-dx))/((x-dx)^3))/(2*dx);
  out = abs(central_differencing_second_order - analytical_derivative);
endfunction

function out = fourth_order_approximation(x,dx)
  %given function f(x) = sin(x)/x^3
  analytical_derivative = (x^3*cos(x)-sin(x)*3*x^2)/x^6;

  %Fourth order approximation using central differencing scheme.
  term1 = (sin(x-(2*dx)))/((x-(2*dx))^3);
  term2 = 8*(sin(x-dx))/((x-dx)^3);
  term3 = 8*(sin(x+dx))/((x+dx)^3);
  term4 = (sin(x+(2*dx)))/((x+(2*dx))^3);
  central_differencing_fourth_order = (term1 - term2 + term3 - term4)/(12*dx);
  out = abs(central_differencing_fourth_order - analytical_derivative);
endfunction

Ouputs:

Result:

• It is clearly evident from the above graph that the convergence of fourth order approximation is quite high i.e. the error reduces in each step.
• Also, we can observe a common trend that with increase in step size (i.e. dx), the absolute error in approximation (irrespective of type) increases.

References:

• Introduction to CFD using MATLAB and OpenFOAM - Course provided by Skill Lync.
• Higher Engineering Mathematics - BS Grewal
• GNU Octave