高斯-勒让德求积工具

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] 纯文本查看 复制代码

查看代码

打印?

001	!--------------------------------------------------------------------------------

002	!EXAMPLE INTEGERAL SIN(X), X FROM 0.0 TO PI

003	!--------------------------------------------------------------------------------

004	PROGRAM TEST_INTEGRAL

005	IMPLICIT NONE

006

007	REAL(KIND=8)::A,B

008	REAL(KIND=8),EXTERNAL::YOUR_FUNC

009	REAL(KIND=8),EXTERNAL::GUALEGINTG12

010	INTEGER,PARAMETER::M=5

011

012	! A1   B1/A2 B2/A3  B3/A4...BM-1/AM  BM

013	! A|-----|-----|------|.......|-----|B

014

015	REAL(KIND=8)::AX(M)

016	REAL(KIND=8)::BX(M)

017	REAL(KIND=8)::DX

018	REAL(KIND=8)::OTHERS=0.D0

019	REAL(KIND=8)::INTGVAL

020	INTEGER::I

021

022	WRITE(*,*)

023	WRITE(*,*)

024	WRITE(*,FMT='(4X,A)',ADVANCE='NO'),'A = '

025	READ(*,FMT='(F15.8)')A

026	WRITE(*,FMT='(4X,A)',ADVANCE='NO'),'B = '

027	READ(*,FMT='(F15.8)')B

028

029	DX=(B-A)/REAL(M,KIND=8)

030	INTGVAL=0.D0

031

DO I=1,M

032	AX(I)=A+REAL(I-1,KIND=8)*DX

033	BX(I)=AX(I)+DX

034	INTGVAL=INTGVAL+GUALEGINTG12(YOUR_FUNC,OTHERS,AX(I),BX(I))

035

ENDDO

036

037	WRITE(*,FMT='(4X,A10,F12.9)'),'INTGVAL = ',INTGVAL

038

READ(*,*)

039

040	END PROGRAM TEST_INTEGRAL

041	!--------------------------------------------------------------------------------

042

043	REAL(KIND=8) FUNCTION YOUR_FUNC(X,OTHERS)

044	IMPLICIT NONE

045	REAL(KIND=8),INTENT(IN)::X,OTHERS

046	!................

047	!FOR EXAMPLE

048	YOUR_FUNC=SIN(X)+OTHERS !OR YOUR_FUNC=X*X+OTHERS

049

RETURN

050	END FUNCTION YOUR_FUNC

051	!--------------------------------------------------------------------------------

052

053	REAL(KIND=8) FUNCTION GUALEGINTG12(YOUR_FUNC,OTHERS,A,B)

054	IMPLICIT NONE

055	REAL(KIND=8),INTENT(IN)::A,B

056	REAL(KIND=8)::YOUR_FUNC,OTHERS

057	EXTERNAL::YOUR_FUNC

058	REAL(KIND=8)::RSLT

059	REAL(KIND=8)::X,APB2,ASB2

060	INTEGER::I

061	!--------------------------------------------------------------------------------

062	!12-POINT GAUSS-LEGENDRE WEIGHTS AND ABSCISSAS

063	!--------------------------------------------------------------------------------

064	!GAUSS WEIGHTS

065	REAL(KIND=8)::WEIGHTS(12)=(/&

066	0.4717533638651182719461596148501741D-01,&

067	0.1069393259953184309602547181939961D+00,&

068	0.1600783285433462263346525295433591D+00,&

069	0.2031674267230659217490644558097984D+00,&

070	0.2334925365383548087608498989248781D+00,&

071	0.2491470458134027850005624360429511D+00,&

072	0.2491470458134027850005624360429511D+00,&

073	0.2334925365383548087608498989248781D+00,&

074	0.2031674267230659217490644558097984D+00,&

075	0.1600783285433462263346525295433591D+00,&

076	0.1069393259953184309602547181939961D+00,&

077	0.4717533638651182719461596148501741D-01 &

078

/)

079	!GAUSS NODES IN INTERVAL[-1.0, 1.0]

080	REAL(KIND=8)::NODES(12)=(/&

081	-0.9815606342467192506905490901492809D+00,&

082	-0.9041172563704748566784658661190962D+00,&

083	-0.7699026741943046870368938332128180D+00,&

084	-0.5873179542866174472967024189405343D+00,&

085	-0.3678314989981801937526915366437176D+00,&

086	-0.1252334085114689154724413694638531D+00,&

087	0.1252334085114689154724413694638531D+00,&

088	0.3678314989981801937526915366437176D+00,&

089	0.5873179542866174472967024189405343D+00,&

090	0.7699026741943046870368938332128180D+00,&

091	0.9041172563704748566784658661190962D+00,&

092	0.9815606342467192506905490901492809D+00 &

093

/)

094	!--------------------------------------------------------------------------------

095

RSLT=0.D0

096	ASB2=0.5D0*(B-A)

097	APB2=0.5D0*(B+A)

098

DO I=1,12

099	X=ASB2*NODES(I)+APB2

100	RSLT=RSLT+WEIGHTS(I)*YOUR_FUNC(X,OTHERS)

101

ENDDO

102	RSLT=RSLT*ASB2

103	GUALEGINTG12=RSLT

104

RETURN

105	END FUNCTION GUALEGINTG12

106	!--------------------------------------------------------------------------------

chiangtp

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

謝謝分享

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

[Fortran] 纯文本查看 复制代码

查看代码

打印?

01	PROGRAM test

02	INTEGER, PARAMETER :: dp=KIND(1.0D0)

03	REAL(KIND=dp) :: x=1.2, y=1.2D0     ! (KIND=8),  literal "8" is compiler dependent

04	WRITE(*,*) x, y                     ! x, default precision only

05	END PROGRAM test

06

07	WRITE(*,'(4X,A)',ADVANCE='NO') 'A = ' ! [), 'A = '] ---> [) 'A = '],  no ","

08	READ(*,*) A                           ! (F15.8), formatted READ, dangerous

09

10	FUNCTION GUALEGINTG12(YOUR_FUNC, OTHERS, A, B)

11	IMPLICIT NONE

12	INTEGER, PARAMETER :: dp=KIND(1.0D0)

13	REAL(KIND=dp)             :: GUALEGINTG12

14	REAL(KIND=dp), EXTERNAL   :: YOUR_FUNC

15	REAL(KIND=dp), INTENT(IN) :: OTHERS, A, B

kerb

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

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

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

ncutefly

不错的资料，谢谢楼主