MATLAB Examples

# Cubic polynomial differintegral

Differintegral of a cubic polynomial using Fourier series representation

## 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.

__________________________________________________________________________
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/

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

__________________________________________________________________________
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 at the range is as follows:

Starting from the following identity:

the order derivative of can be written as:

From this result, it is clear that when order goes from to , the derivative goes from to by gradually increasing amplitude from to and phase from to . 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 order derivative of the function defined above is calculated as:

where the constant term 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:2*pi)';
%2
f=@(x) x.*(x-pi).*(x-2*pi);
%3
a=x(1);
b=x(end);
%4
numFourier=60;
%5
numGauss=120;
%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)
4. Calculate the first integral of the function using its Fourier series representation (diffintF-1)
%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
diffintF05=fourier_diffint(f,x,deg1,numFourier,numGauss);
%4
diffintF_1=fourier_diffint(f,x,deg2,numFourier,numGauss);

## Plots

1. Plot the function using its definition and its Fourier series representation (check accuracy)
2. Plot the half-derivative of the function using its Fourier series representation (diffintF0.5)
3. Plot the first integral of the function using its Fourier series representation (diffintF-1)
%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,diffintF05,'LineWidth',2.)
grid on
xlabel('x','FontSize',13);
ylabel('F _ 0.5','FontSize',13);
title('Half-derivative','FontSize',13)
%3
figure('Name','First integral','NumberTitle','off')
plot(x,diffintF_1,'LineWidth',2.)
grid on
xlabel('x','FontSize',13);
ylabel('F _ -1','FontSize',13);
title('First integral','FontSize',13)