MATLAB Examples

Identity function differintegral

Differintegral of the identity function using Fourier series representation

Contents

Notation and references

The notation followed here and in the following MATLAB codes:

  • fourier_diffint.m
  • fourier.m

does not conform to a specific reference, but it is explained into the help sections of the above codes.

Functions used

  • Calculation of Fourier series coefficients
  • Calculation of function differintegral based on Fourier series representation
help fourier
help fourier_diffint
  Fourier series of an arbitrary function.
 
  Description
      The coefficients of the Fourier series are extracted, including a0,
      for a function defined in a given range [#a# #b#]. The necessary
      integrations are performed with the Gauss-Legendre quadrature rule.
      Selection for the number of desired Fourier coefficient pairs is made
      as well as for the number of the Gauss-Legendre integration points.
      Unlike many publicly available functions, this function can work for
      #numGauss#>=46. It does not rely on the build-in Matlab routine
      'roots' to determine the roots of the Legendre polynomial, but finds
      the roots by looking for the eigenvalues of an alternative version of
      the companion matrix of the n'th degree Legendre polynomial. The
      companion matrix is constructed as a symmetrical matrix, guaranteeing
      that all the eigenvalues (roots) will be real. On the contrary, the
      'roots' function uses a general form for the companion matrix, which
      becomes unstable at higher values of #numGauss#, leading to complex
      roots.
 
  Required input arguments
      #f# (function handle) defines the function to be developed into
      Fourier series. The function must use array operators instead of
      matrix operators. For example, use .* (times) rather than * (mtimes).
      #a# (scalar) is the lower limit of the range in which #f# is defined.
      #b# (scalar) is the upper limit of the range in which #f# is defined.
      #numFourier# (scalar) is the number of Fourier coefficients desired,
      including a0.
      #numGauss# (scalar) is the number of Gauss integration points used.
      Experience has shown that there must be #numGauss#>=2*#numFourier# in
      order to achieve adequately accurate results.
 
  Output parameters
      #an# ([#numFourier# x 1]) contains the alpha coefficients of the
      Fourier series.
      #bn# ([#numFourier# x 1]) contains the beta coefficients of the
      Fourier series.
 
  Parents (calling functions)
      fourier_diffint.m
 
  Children (called functions)
      None.
 
 __________________________________________________________________________
  Copyright (c) 13-Mar-2014
      George Papazafeiropoulos
      First Lieutenant, Infrastructure Engineer, Hellenic Air Force
      Civil Engineer, M.Sc., Ph.D. candidate, NTUA
      Email: gpapazafeiropoulos@yahoo.gr
      Website: http://users.ntua.gr/gpapazaf/
 
 

    Overloaded methods:
       sym/fourier

  Fractional differintegral of an arbitrary function.
 
  Description
      The #deg#-th order differintegral of the function #f# defined in a
      given range [#a# #b#] is calculated through Fourier series expansion,
      where #deg# is any real number and not necessarily integer. The
      necessary integrations are performed with the Gauss-Legendre
      quadrature rule. Selection for the number of desired Fourier
      coefficient pairs is made as well as for the number of the
      Gauss-Legendre integration points.
 
  Required input arguments
      #f# (function handle) defines the function to be developed into
      Fourier series. The function must use array operators instead of
      matrix operators. For example, use .* (times) rather than * (mtimes).
      #x# ([#rangelength# x 1]) is the vector defining the differentiation
      or integration range.
      #deg# (scalar) is the order of differentiation.
      #numFourier# (scalar) is the number of Fourier coefficients desired,
      including a0.
      #numGauss# (scalar) is the number of Gauss integration points used.
      Experience has shown that there must be #numGauss#>=2*#numFourier# in
      order to achieve adequately accurate results.
 
  Output parameters
      #diffintF# ([#rangelength# x 1]) is the fractional differintegral at
      the range #x#.
 
  Parents (calling functions)
      None.
 
  Children (called functions)
      fourier.m
 
 __________________________________________________________________________
  Copyright (c) 13-Mar-2014
      George Papazafeiropoulos
      First Lieutenant, Infrastructure Engineer, Hellenic Air Force
      Civil Engineer, M.Sc., Ph.D. candidate, NTUA
      Email: gpapazafeiropoulos@yahoo.gr
      Website: http://users.ntua.gr/gpapazaf/
 
 

Theory

The Fourier series representation of a function $$F$ at the range $$\left[a,b\right]$ is as follows:

$$F=\mathrm{A_0} + \sum_{k=1}^{\infty} \mathrm{A_k}\, \cos\left(\frac{2\,
\pi\, k\, x}{b - a}\right) + \sum_{k=1}^{\infty} \mathrm{B_k}\,
\sin\left(\frac{2\, \pi\, k\, x}{b - a}\right)$$

Starting from the following identity:

$$\frac{d\left[sin(a\, x)\right]}{dx}=a\, cos(a\, x)=a\, sin(a\,
x+ \frac{\pi}{2})$$

the $$n^{th}$ order derivative of $$sin(a\, x)$ can be written as:

$$\frac{d^{n}\left[sin(a\, x)\right]}{dx^{n}}=a^{n}\, sin(a\, x+
\frac{n\, \pi}{2})$$

From this result, it is clear that when order $$n$ goes from $$0$ to $$1$, the derivative $$\frac{d^{n}\left[sin(a\, x)\right]}{dx^{n}}$ goes from $$sin(a\, x)$ to $$a\, cos(a\, x)=a\, sin(a\, x+\pi/2)$ by gradually increasing amplitude from $$1$ to $$a$ and phase from $$0$ to $$\pi/2$. Thus, the fractional derivative can be interpreted as a function interpolation between derivatives with integer order.

The same identities as the above hold for the cosine function. Based on this rationale, the $$n^{th}$ order derivative of the function $$F$ defined above is calculated as:

$$\frac{d^{n}F}{dx^{n}}=\mathrm{A_0}\, \frac{x^{-n}}{\Gamma(1-n)} +
\sum_{k=1}^{\infty} \mathrm{A_k}\, \left(\frac{2\, \pi\, k}{b -
a}\right)^{n}\, \cos\left(\frac{2\, \pi\, k\, x}{b - a}+\frac{n\, \pi}{2}
\right) + \sum_{k=1}^{\infty} \mathrm{B_k}\, \left(\frac{2\, \pi\, k}{b -
a}\right)^{n}\, \sin\left(\frac{2\, \pi\, k\, x}{b - a}+\frac{n\, \pi}{2}
\right)$$

where the constant term $$A_0$ of the Fourier series is differentiated according to well-known relations of fractional calculus.

Initial definitions

In the subsequent code the following initial definitions are made (in the order presented below):

  1. Set the range of differintegration
  2. Define the function to be differintegrated
  3. Set the range limits
  4. Set the number of Fourier terms
  5. Set the number of Gauss points
  6. Set the degree of half-derivative
  7. Set the degree of first integral
%1
x=(0:0.01:1)';
%2
f=@(x) x;
%3
a=x(1);
b=x(end);
%4
numFourier=260;
%5
numGauss=520;
%6
deg1=0.5;
%7
deg2=-1;

Applications

  1. Calculate the function using its definition
  2. Calculate the function using its Fourier series representation and check accuracy
  3. Calculate the half-derivative of the function using its Fourier series representation (diffintF0.5,1)
  4. Calculate the half-derivative of the function using its definition (diffintF0.5,2)
  5. Calculate the first integral of the function using its Fourier series representation (diffintF-1). It is equal to the area under function f
%1
F1=f(x);
%2
[A,B]=fourier(f,a,b,numFourier,numGauss);
n=(0:numFourier-1)';
n=n(:,ones(1,length(x)));
An=A(:,ones(1,length(x)));
Bn=B(:,ones(1,length(x)));
evalx=x(:,ones(1,numFourier))';
F2=sum(An.*cos(2*pi*n.*evalx/(b-a))+Bn.*sin(2*pi*n.*evalx/(b-a)))';
%3
diffintF051=fourier_diffint(f,x,deg1,numFourier,numGauss);
%4
diffintF052=@(x) 2*sqrt(x/pi);
%5
diffintF_1=fourier_diffint(f,x,deg2,numFourier,numGauss);
area=diffintF_1(end)-diffintF_1(1);
zero=area-1^2/2
zero =

  -1.2212e-15

Plots

  1. Plot the function using its definition and its Fourier series representation and check accuracy
  2. Plot the half-derivative of the function using its Fourier series representation and its definition (diffintF0.5) and check accuracy
%1
figure('Name','Function plot','NumberTitle','off')
plot(x,F1,'LineWidth',2.)
grid on
xlabel('x','FontSize',13);
ylabel('F1,F2','FontSize',13);
title('Function plot','FontSize',13)
hold on
plot(x,F2,'LineWidth',2.)
hold off
%2
figure('Name','Half-derivative','NumberTitle','off')
plot(x,diffintF051,'LineWidth',2.)
grid on
xlabel('x','FontSize',13);
ylabel('F _ 0.5_1,F _ 0.5_2','FontSize',13);
title('Half-derivative','FontSize',13)
hold on
plot(x,diffintF052(x),'LineWidth',2.)
hold off

Copyright

Copyright (c) 13-Mar-2014 by George Papazafeiropoulos