MATLAB Examples

# Identity function differintegral

Differintegral of the identity function 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. __________________________________________________________________________ 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 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: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