高斯-勒让德求积工具

kerb

看见论坛上有相关的帖子，觉得那样并不方便，于是就做了个工具，高斯-勒让德权重系数和节点计算的代码代码来自numerica recipes in fortran77

获得较高积分精度并不是简单增加阶数就可以的，好的做法是把积分区间分成适当的份数，选用适当阶数的高斯节点

高斯-勒让德求积工具
输入：勒让德多项式的阶数
输出：
         屏幕是计算好的权重系数和勒让德节点
         同时文件输出一个fortran代码的例子，比如输入4，文件内容如下
FORTRAN版高斯乐让德积分工具 GAULEGQUADF.zip (275.85 KB, 下载次数: 18)

2017-9-21 18:18 上传

点击文件名下载附件

!--------------------------------------------------------------------------------
!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   
!--------------------------------------------------------------------------------[/mw_shl_code]

kerb

第二版
第二版会生成一个例子程序，例子是sin(x)在区间[0.0,3.1415926]的积分，区间分成5段（程序中整形变量M）
FORTRAN版高斯乐让德积分工具 GAULEGQUADF.zip (275.85 KB, 下载次数: 15)

2017-9-21 18:19 上传

点击文件名下载附件

这是输入12后生成的

以下是例子代码

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

!--------------------------------------------------------------------------------
!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/A3  B3/A4...BM-1/AM  BM
    ! 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   
!--------------------------------------------------------------------------------

chiangtp

kerb 发表于 2017-9-17 21:38
第二版
第二版会生成一个例子程序，例子是sin(x)在区间[0.0,3.1415926]的积分，区间分成5段（程序中整形变 ...

謝謝分享

雞蛋裡挑個骨頭, 還請包含

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

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

chiangtp 发表于 2017-9-18 13:32
謝謝分享

雞蛋裡挑個骨頭, 還請包含

谢谢提醒，等我有时间在进一步完善
已改正

ncutefly

不错的资料，谢谢楼主