高斯-勒让德求积工具
本帖最后由 kerb 于 2017-9-21 18:18 编辑看见论坛上有相关的帖子，觉得那样并不方便，于是就做了个工具，高斯-勒让德权重系数和节点计算的代码代码来自numerica recipes in fortran77
获得较高积分精度并不是简单增加阶数就可以的，好的做法是把积分区间分成适当的份数，选用适当阶数的高斯节点
高斯-勒让德求积工具
输入：勒让德多项式的阶数
输出：
屏幕是计算好的权重系数和勒让德节点
同时文件输出一个fortran代码的例子，比如输入4，文件内容如下
!--------------------------------------------------------------------------------
!EXAMPLE
!--------------------------------------------------------------------------------
REAL(KIND=8) FUNCTION YOUR_FUNC(X,OTHERS)
IMPLICIT NONE
REAL(KIND=8),INTENT(IN)::X,OTHERS
!................
!FOR EXAMPLE
YOUR_FUNC=X*X+OTHERS
RETURN
END FUNCTION YOUR_FUNC
!--------------------------------------------------------------------------------
REAL(KIND=8) FUNCTION GUALEGINTG4 (YOUR_FUNC,OTHERS,A,B)
IMPLICIT NONE
REAL(KIND=8),INTENT(IN)::A,B
REAL(KIND=8)::YOUR_FUNC,OTHERS
EXTERNAL::YOUR_FUNC
REAL(KIND=8)::RSLT
REAL(KIND=8)::X,APB2,ASB2
INTEGER::I
!4 -POINT GAUSS-LEGENDRE WEIGHTS AND ABSCISSAS
!--------------------------------------------------------------------------------
!GAUSS WEIGHTS
REAL(KIND=8)::WEIGHTS(4 )=(/&
0.3478548451374538573730639492219996 ,&
0.6521451548625461426269360507780006 ,&
0.6521451548625461426269360507780006 ,&
0.3478548451374538573730639492219996 &
/)
!GAUSS NODES IN INTERVAL[-1.0, 1.0]
REAL(KIND=8)::NODES(4 )=(/&
-0.8611363115940525752239464888928096 ,&
-0.3399810435848562648026657591032447 ,&
0.3399810435848562648026657591032447 ,&
0.8611363115940525752239464888928096 &
/)
!--------------------------------------------------------------------------------
RSLT=0.D0
ASB2=0.5D0*(B-A)
APB2=0.5D0*(B+A)
DO I=1,4
X=ASB2*NODES(I)+APB2
RSLT=RSLT+WEIGHTS(I)*YOUR_FUNC(X,OTHERS)
ENDDO
RSLT=RSLT*ASB2
GUALEGINTG4 =RSLT
RETURN
END FUNCTION GUALEGINTG4
!--------------------------------------------------------------------------------
本帖最后由 kerb 于 2017-9-21 18:19 编辑
第二版
第二版会生成一个例子程序，例子是sin(x)在区间的积分，区间分成5段（程序中整形变量M）
这是输入12后生成的
以下是例子代码
!--------------------------------------------------------------------------------
!EXAMPLE INTEGERAL SIN(X), X FROM 0.0 TO PI
!--------------------------------------------------------------------------------
PROGRAM TEST_INTEGRAL
IMPLICIT NONE
REAL(KIND=8)::A,B
REAL(KIND=8),EXTERNAL::YOUR_FUNC
REAL(KIND=8),EXTERNAL::GUALEGINTG12
INTEGER,PARAMETER::M=5
! A1 B1/A2 B2/A3B3/A4...BM-1/AMBM
! A|-----|-----|------|.......|-----|B
REAL(KIND=8)::AX(M)
REAL(KIND=8)::BX(M)
REAL(KIND=8)::DX
REAL(KIND=8)::OTHERS=0.D0
REAL(KIND=8)::INTGVAL
INTEGER::I
WRITE(*,*)
WRITE(*,*)
WRITE(*,FMT='(4X,A)',ADVANCE='NO'),'A = '
READ(*,FMT='(F15.8)')A
WRITE(*,FMT='(4X,A)',ADVANCE='NO'),'B = '
READ(*,FMT='(F15.8)')B
DX=(B-A)/REAL(M,KIND=8)
INTGVAL=0.D0
DO I=1,M
AX(I)=A+REAL(I-1,KIND=8)*DX
BX(I)=AX(I)+DX
INTGVAL=INTGVAL+GUALEGINTG12(YOUR_FUNC,OTHERS,AX(I),BX(I))
ENDDO
WRITE(*,FMT='(4X,A10,F12.9)'),'INTGVAL = ',INTGVAL
READ(*,*)
END PROGRAM TEST_INTEGRAL
!--------------------------------------------------------------------------------
REAL(KIND=8) FUNCTION YOUR_FUNC(X,OTHERS)
IMPLICIT NONE
REAL(KIND=8),INTENT(IN)::X,OTHERS
!................
!FOR EXAMPLE
YOUR_FUNC=SIN(X)+OTHERS !OR YOUR_FUNC=X*X+OTHERS
RETURN
END FUNCTION YOUR_FUNC
!--------------------------------------------------------------------------------
REAL(KIND=8) FUNCTION GUALEGINTG12(YOUR_FUNC,OTHERS,A,B)
IMPLICIT NONE
REAL(KIND=8),INTENT(IN)::A,B
REAL(KIND=8)::YOUR_FUNC,OTHERS
EXTERNAL::YOUR_FUNC
REAL(KIND=8)::RSLT
REAL(KIND=8)::X,APB2,ASB2
INTEGER::I
!--------------------------------------------------------------------------------
!12-POINT GAUSS-LEGENDRE WEIGHTS AND ABSCISSAS
!--------------------------------------------------------------------------------
!GAUSS WEIGHTS
REAL(KIND=8)::WEIGHTS(12)=(/&
0.4717533638651182719461596148501741D-01,&
0.1069393259953184309602547181939961D+00,&
0.1600783285433462263346525295433591D+00,&
0.2031674267230659217490644558097984D+00,&
0.2334925365383548087608498989248781D+00,&
0.2491470458134027850005624360429511D+00,&
0.2491470458134027850005624360429511D+00,&
0.2334925365383548087608498989248781D+00,&
0.2031674267230659217490644558097984D+00,&
0.1600783285433462263346525295433591D+00,&
0.1069393259953184309602547181939961D+00,&
0.4717533638651182719461596148501741D-01 &
/)
!GAUSS NODES IN INTERVAL[-1.0, 1.0]
REAL(KIND=8)::NODES(12)=(/&
-0.9815606342467192506905490901492809D+00,&
-0.9041172563704748566784658661190962D+00,&
-0.7699026741943046870368938332128180D+00,&
-0.5873179542866174472967024189405343D+00,&
-0.3678314989981801937526915366437176D+00,&
-0.1252334085114689154724413694638531D+00,&
0.1252334085114689154724413694638531D+00,&
0.3678314989981801937526915366437176D+00,&
0.5873179542866174472967024189405343D+00,&
0.7699026741943046870368938332128180D+00,&
0.9041172563704748566784658661190962D+00,&
0.9815606342467192506905490901492809D+00 &
/)
!--------------------------------------------------------------------------------
RSLT=0.D0
ASB2=0.5D0*(B-A)
APB2=0.5D0*(B+A)
DO I=1,12
X=ASB2*NODES(I)+APB2
RSLT=RSLT+WEIGHTS(I)*YOUR_FUNC(X,OTHERS)
ENDDO
RSLT=RSLT*ASB2
GUALEGINTG12=RSLT
RETURN
END FUNCTION GUALEGINTG12
!--------------------------------------------------------------------------------
kerb 发表于 2017-9-17 21:38
第二版
第二版会生成一个例子程序，例子是sin(x)在区间的积分，区间分成5段（程序中整形变 ...
謝謝分享
雞蛋裡挑個骨頭, 還請包含
PROGRAM test
INTEGER, PARAMETER :: dp=KIND(1.0D0)
REAL(KIND=dp) :: x=1.2, y=1.2D0 ! (KIND=8),literal "8" is compiler dependent
WRITE(*,*) x, y ! x, default precision only
END PROGRAM test
WRITE(*,'(4X,A)',ADVANCE='NO') 'A = ' ! [), 'A = '] ---> [) 'A = '],no ","
READ(*,*) A ! (F15.8), formatted READ, dangerous
FUNCTION GUALEGINTG12(YOUR_FUNC, OTHERS, A, B)
IMPLICIT NONE
INTEGER, PARAMETER :: dp=KIND(1.0D0)
REAL(KIND=dp) :: GUALEGINTG12
REAL(KIND=dp), EXTERNAL :: YOUR_FUNC
REAL(KIND=dp), INTENT(IN) :: OTHERS, A, B 本帖最后由 kerb 于 2017-9-21 18:20 编辑
chiangtp 发表于 2017-9-18 13:32
謝謝分享
雞蛋裡挑個骨頭, 還請包含
谢谢提醒，等我有时间在进一步完善
已改正
页:
[1]