WANtaroHP (F90 MATH)

Outline of this page

Contents
Simultaneous Linear Equations
Eigenvalue Analysis
Gauss-Legendre rule
Fast Fourier Transform (FFT and Inverse FFT)
Sample use of Complex number
Bessel functions, Gamma function (Fortran 2008)
Cubic Spline Interpolation
Solution of Algebraic Equations (DKA method)
Complete elliptic integral

Simultaneous Linear Equations

Outline of this program

This is a program to solve simultaneous linear equations.
Gauss-Jordan elimination and Cholesky method for band matrix are included in this program. These programs are used in FEM programs as a subroutine.

Source code by f90

Filename	Description
f90_GJ.txt	program to solve simultaneous linear equations

Eigenvalue Analysis

Outline of this program

This is a program for eigenvalue analysis.
Jacobi method and generalized Jacobi method are included in this program.
Jacobi method is used in principal component analysis as a subroutine.
Generalized Jacobi method is used in FEM programs for buckling analysis and eigenfrequency analysis as a subroutine.

Source code by f90

Filename	Description
f90_JACOBI.txt	program for eigenvalue analysis

Gauss-Legendre rule

Outline of this program

In this section, two programs are introduced.
One is a program to calculate the coordinates and weights of Gauss point.
Another one is a program which shows the adoption sample of Gauss-Legendre rule.
The document of 'http://www7.ocn.ne.jp/~kawa1/numeric.pdf' was referred to make my programs.

Source code by f90

Filename	Description
f90_GAUSSIP.txt	program to calculate the coordinates and weights of Gauss point
f90_CALPI.txt	program which shows the adoption sample of Gauss-Legendre rule

Fast Fourier Transform (FFT and Inverse FFT)

Outline of this program

This is a program for FFT and IFFT.
The complex variables on Fortran are not used in this program. Calculation is carried out by preparing the array for real parts and the array for imaginary parts separately.
This program is used in the programs for Fourier spectrum analysis etc. as a subroutine.

Source code by f90

Filename	Description
f90_FFT.txt	program for FFT and IFFT
inp_FFT_0.txt	Sample input data
out_FFT_0.txt	Sample output data

Sample use of Complex number

Outline of this program

This is a program for sample use of complex number.

Source code by f90

Filename	Description
f90_COMPLEX.txt	Program for sample use of complex number

Bessel functions, Gamma function (Fortran 2008)

Outline of this program

This is a program for sample use of Bessel functions and Gamma function in Fortran 2008.
Png image files by gnuplot are also shown. In png image files, the plots using gnuplot functions are included for references.
Bach command for execution are shown below.

gfortran -o f08_BESSEL.exe f08_BESSEL.f08
gfortran -o f08_GAMMA.exe f08_GAMMA.f08
gfortran -o f08_GPL_PLOT.exe f08_GPL_PLOT.f08

f08_BESSEL > test_b.txt
f08_GAMMA > test_g.txt
f08_GPL_PLOT test_b.txt test_g.txt

Source code by f08

Filename	Description
f08_BESSEL.txt	Program for sample use of Bessel functions
f08_GAMMA.txt	Program for sample use of Gamma function
f08_GPL_PLOT.txt	Program for gnuplot
fig_png_bessel.png	png image file of Bessel functions plot
fig_png_gamma.png	png image file of Gamma functions plot

Cubic Spline Interpolation

Outline of this program

This is a program for Cubic Spline Interpolation.
As a solver for simultaneous equations, Cholesky method for band matrix is adopted.
Script for gnuplot drawing is also attached in this section because of quick confirmation of the result.

Bach command for execution

gfortran -o f90_SPL.exe f90_SPL.f90
f90_SPL inp_data.txt 5 > out_data.txt

The format of command line arguments are shown below.

f90_SPL fnameR md > fnameW

f90_SPL	Compiled F90 frogram for Cubic Spline Interpolation
fnameR	Input file name
md	Number of division between the one plot to next plot
fnameW	Output file name

Source code by f90

Filename	Description
f90_SPL.txt	Program for Cubic Spline Interpolation
inp_data.txt	Input data sample
out_data.txt	Output data sample
gpl_spline.txt	Script for gnuplot drawing
fig_png_spline.png	png image file of the result

Solution of Algebraic Equations (DKA method)

Outline of this program

This is a program for solution of algebraic equations. This can also obtain complex number solutions.
Solution method is DKA method. This was developed by Durand, Kerner and Aberth.
The quadruple precision type is used for the calculation, and the number of significant digits is about 30.
When multiple root is included in the solution, solution becomes low accuracy.
In this program, complex type variants are not used, and real part and imaginary part are treated separately.
This Fortran program was made refering to the C program shown in following document.

http://www.ecs.shimane-u.ac.jp/~kyoshida/c6(2002).pdf

Bach command for execution

gfortran -o f90_DKA.exe f90_DKA.f90

f90_DKA inp_test10.txt > out_test10.txt
copy _data_circ.txt _data_circ10.txt
copy _gmt_inp.txt _gmt_inp10.txt

f90_DKA inp_test20.txt > out_test20.txt
copy _data_circ.txt _data_circ20.txt
copy _gmt_inp.txt _gmt_inp20.txt

call bat_gmt_draw

The format of command line arguments are shown below.

f90_DKA fnameR > fnameW

f90_DKA	Compiled F90 program
fnameR	Input file name
fnameW	Output file name

The purpose of output file 'fnameW' is to confirm the numerical results. In addition, following files which have fixed file name for each are created for GMT drawing.

_iteration.txt	Work file
_gmt_inp.txt	Data file for GMT (coordinate data of coefficients included convergence process)
_data_circ.txt	Data file for GMT ( coordinate data for a circle related Aberth's initial values)

Source code by f90

Filename	Description
f90_DKA.txt	Program for DKA method
inp_test10.txt	Input data sample (Case-1)
inp_test20.txt	Input data sample (Case-2)
inp_test30.txt	Input data sample (Case-3)
out_test10.txt	Output data sample (Case-1)
out_test20.txt	Output data sample (Case-2)
out_test30.txt	Output data sample (Case-3)
bat_gmt_draw.txt	Batch filr for GMT drawing
fig_gmt_out10_1.png	png image file (Case-1)
fig_gmt_out10_2.png	png image file (Case-1 detail)
fig_gmt_out20_1.png	png image file (Case-2)
fig_gmt_out20_2.png	png image file (Case-2 detail)
fig_gmt_out30_1.png	png image file (Case-3)
fig_gmt_out30_2.png	png image file (Case-3 detail)
TeX_DKA.pdf	Document

Case-1	f(x)=(x+1)^10=0
Case-2	f(x)=(x-20)(x-19)...(x-1)=0
Case-3	f(x)=(x-30)(x-29)...(x-1)=0

Complete elliptic integral

Outline of this program

These are programs for Complete elliptic integrals of the 1st kind K(p) and 2nd kind E(p).
The numerical integration methods are one is by Series expansion and the other is by Gauss-Legendre rule.
The method by Gauss-Legendre rule is faster than that by Series expansion.
Samples of calculation results are shown below. ita. means number for conversion for each calculation.

Result by Series expansion method

  Mehod by Series expansion
  ita.   p              K(p)           E(p)
    1  0.0000000E+00  0.1570796E+01  0.1570796E+01
    5  0.1000000E+00  0.1574746E+01  0.1566862E+01
    7  0.2000000E+00  0.1586868E+01  0.1554969E+01
    9  0.3000000E+00  0.1608049E+01  0.1534833E+01
   11  0.4000000E+00  0.1640000E+01  0.1505942E+01
   14  0.5000000E+00  0.1685750E+01  0.1467462E+01
   19  0.6000000E+00  0.1750754E+01  0.1418083E+01
   27  0.7000000E+00  0.1845694E+01  0.1355661E+01
   41  0.8000000E+00  0.1995303E+01  0.1276350E+01
   83  0.9000000E+00  0.2280549E+01  0.1171697E+01
54842  0.9999000E+00  0.5645147E+01  0.1000515E+01
 time=   47.875000     (sec)

Result by Gauss-Legendre rule: number of integration = 10

  method by Gauss-Legendre rule (n=10)
  ita.   p              K(p)           E(p)
    2  0.0000000E+00  0.1570796E+01  0.1570796E+01
    2  0.1000000E+00  0.1574746E+01  0.1566862E+01
    2  0.2000000E+00  0.1586868E+01  0.1554969E+01
    2  0.3000000E+00  0.1608049E+01  0.1534833E+01
    2  0.4000000E+00  0.1640000E+01  0.1505942E+01
    2  0.5000000E+00  0.1685750E+01  0.1467462E+01
    2  0.6000000E+00  0.1750754E+01  0.1418083E+01
    2  0.7000000E+00  0.1845694E+01  0.1355661E+01
    3  0.8000000E+00  0.1995303E+01  0.1276350E+01
    3  0.9000000E+00  0.2280549E+01  0.1171697E+01
   53  0.9999000E+00  0.5645148E+01  0.1000515E+01
 time=   0.0000000     (sec)

Result by Gauss-Legendre rule: number of integration = 2

  method by Gauss-Legendre rule (n=2)
  ita.   p              K(p)           E(p)
    2  0.0000000E+00  0.1570796E+01  0.1570796E+01
    3  0.1000000E+00  0.1574746E+01  0.1566862E+01
    4  0.2000000E+00  0.1586868E+01  0.1554969E+01
    4  0.3000000E+00  0.1608049E+01  0.1534833E+01
    5  0.4000000E+00  0.1640000E+01  0.1505942E+01
    6  0.5000000E+00  0.1685750E+01  0.1467462E+01
    6  0.6000000E+00  0.1750754E+01  0.1418083E+01
    7  0.7000000E+00  0.1845694E+01  0.1355661E+01
    9  0.8000000E+00  0.1995303E+01  0.1276350E+01
   13  0.9000000E+00  0.2280549E+01  0.1171697E+01
  323  0.9999000E+00  0.5645148E+01  0.1000515E+01
 time=  0.12500000     (sec)

Source code by f90

Filename	Description
f90_CEIS.txt	Program by Series expansion
f90_CEIG.txt	Program by Gauss-Legendre rule