高斯勒让德求积分Fortran程序—自带求任意阶高斯点及权重

岸边的鱼

高斯勒让德求积分Fortran程序——自带求任意阶高斯点及相应权重的程序

句号大师（王培杰）编写，岸边的鱼（小超超小鱼）整理，标注

最近，因为要求解一些剧烈变化（且积分区间跨度比较大）的函数积分问题，想到了高斯型求积分公式，随查阅各种算法书籍，但无奈算法书上给出的都是利用五阶高斯勒让德公式求解积分的程序。经石子儿兄弟指点，认识到五阶高斯勒让德公式求解精度实在是有点小低，（解决办法有两个，一是本文的增加高斯结点，而是采用变步长高斯求积分法，此法将另文讨论。）石子儿兄弟建议将精度至少提高到300—500之间，那将面对自己计算高斯点及权系数的问题。无奈之下，遍寻各大编程网站，终于在我们亲爱的fcode.cn网站上找到了句号大师的一篇帖子，随私聊句号大师,感谢慷慨的句号大师将修改好的程序无私的传给我，感动之余，觉得有必要将这一成果整理，标注后发回我们的网站上供广大网友下载使用。

声明：下文程序整体由句号大师（王培杰）传递给我，仅对其表示由衷的感谢。我负责添加标注，设置了两个输出到文件的功能。另，句号大师提醒广大网友，高斯型求积公式不是高斯点越多越精确，高斯点变多后会有很多的截断误差。

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

module GSLD     !高斯—勒让德积分高斯点及权重的求解模块
implicit none 
    integer,parameter::n=5  !设置求解高斯点的个数
integer,parameter::nn=16 !设置kind的数值
contains
real(kind=nn) function f(x)  !定义四阶精度的函数f（x）
implicit none
  integer::i
  real(kind=nn)::a(n),x      !a(n)代表n阶勒让德多项式
a(1)=x                     !1阶勒让德多项式
a(2)=1.5_16*(x**2)-0.5_16  !2阶勒让德多项式
do i=3,n            
       a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i    
　　　　　　　　　　　　!利用递推关系产生n阶勒让德多项式
enddo
       f=a(n)   !生成的n阶勒让德多项式   
end  function

real(kind=nn) function f1(x)     
  implicit  none
  integer::i
  real(kind=nn)::a(n),x
    a(1)=x
    a(2)=1.5_16*x**2-0.5_16
    do i=3,n-1
       a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
    enddo
       f1=a(n-1)  !生成的（n-1）阶勒让德多项式     
endfunction

real(kind=nn) function g(x)
  implicit none
  integer::i
  real(kind=nn)::a(n),x
    a(1)=x
    a(2)=1.5_16*x**2-0.5_16
    do i=3,n
       a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
    enddo
       g=n*a(n-1)/(1-x**2)-n*x*a(n)/(1-x**2)   
　　　　　　　　　　　　!生成n阶勒让德多项式的导数表达式
endfunction

real(kind=nn) function bis(a,b)   !二分法求解函数的解
implicit  none
  real(kind=nn)::a,b,c       
　　　　　　　　!a,b是传递进来的划分好的有一个解存在的区间
do while(.true.)
      c=(a+b)/2.0_16
      if(f(c)*f(a)<0)then
        b=c
        else
        a=c
      endif
       if((b-a)<1E-34)exit    !这个精度设置的。。（有点狠了）
enddo
       bis=c                  !bis即是利用二分法求得的解
end   function

subroutine fn0(fn,ak)      
  implicit  none
  real(kind=nn)::m,fn(n),ak(n)  
　　　　　　　　　　　　　　!定义数组,大小n由module开始声明。
integer::i,j
open(10,file='高斯点+权重.dat')
    j=0                 !赋值控制循环变量的初值           
m=-1.000001         !设置计算域[-1，1] 的下限，即代替-1 
    do i=1,200000       !这个循环次数应该是由步长0.00001决 定,计算方法：200000=2/0.000001     
      if(f(m)*f(m+0.00001)<0)then !从下限处开始往上逐步累加，
!由步长0.00001说明最多求解10^5个解
j=j+1                         !记录这是第几个解
fn(j)=bis(m,m+0.00001)        
!调用二分法求解程序在分好的一小段上求解，
!将解存储在fn（j）
ak(j)=2.0_16/(n*f1(fn(j))*g(fn(j)))  !高斯点的权重
write(*,*)   '高斯点序号',j
         write(*,*)   '高斯点'    ,fn(j)
         write(*,*)   '高斯点权重',ak(j)
          write(10,*) '高斯点',fn(j)
          write(10,*) '权重'  ,ak(j)
          endif
         m=m+0.00001       !执行完一次判断m向前推进一步
enddo
close (10)
endsubroutine

real(kind=nn) function y(x)  !被积函数
implicit  none
  real(kind=nn)::x
    y=x**2/3.0_16      !每次计算在这里修改被积函数即可
end  function      
end  module

program www_fcode_cn
  use GSLD
  implicit none
  real(kind=nn)::fn(n),ak(n),x,a,b,answer
  integer::i
  open(11,file='GS_LD_result.dat')
  call fn0(fn,ak)      !调用求高斯零点和权重的子函数
a=0.0_16          !积分上限
b=1.0_16          !积分下限   
answer=0.0_16     !求积分结果赋初始值
do i=1,n
        answer=answer+ak(i)*y((a+b)/2.0_16+(b-a)/2.0_16*fn(i)) 
　　　　　　　　　　　　　　  !部分高斯勒让德求积分公式
enddo
        answer=answer*(b-a)/2.0_16
        write(11,*)   '高斯_勒让德求积分结果',answer
        write(*,*)    '高斯_勒让德求积分结果',answer
  close(11)
     stop
end  program www_fcode_cn

运行结果：

以下是使用高斯勒让德函数计算三重积分的例子：

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

module GSLD
    implicit none
    integer,parameter::n=5
    contains

    real*16 function f(x)
    integer::i
    real*16::a(n),x
    a(1)=x
    a(2)=1.5*(x**2)-0.5
    do i=3,n
        a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
    enddo
    f=a(n)
    endfunction

    real*16 function f1(x)
    integer::i
    real*16::a(n),x
    a(1)=x
    a(2)=1.5*x**2-0.5
    do i=3,n
        a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
    enddo
    f1=a(n-1)
    endfunction

    real*16 function g(x)
    integer::i
    real*16::a(n),x
    a(1)=x
    a(2)=1.5d0*x**2-0.5d0
    do i=3,n
        a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
    enddo
    g=n*a(n-1)/(1-x**2)-n*x*a(n)/(1-x**2)
    endfunction

    real*16 function bis(a,b)
    real*16::a,b,c
    do 
        c=(a+b)/2d0
        if(f(c)*f(a)<0)then
            b=c
        else
            a=c
        endif
        if((b-a)<1E-34)exit
    enddo
    bis=c
    endfunction

    subroutine fn0(fn,ak)
    real*16::m,fn(n),ak(n)
    integer::i,j
    j=0
    m=-1.00000001
    do i=1,200000
        if(f(m)*f(m+0.00001d0)<0)then
            j=j+1
            fn(j)=bis(m,m+0.00001)
            ak(j)=2/(n*f1(fn(j))*g(fn(j)))
            write(*,*)j,fn(j),ak(j)
        endif
        m=m+0.00001
    enddo
    endsubroutine

    real*16 function y(x,x1,x2)  !被积函数
real*16::x,x1,x2
    y=2*x+2*x1+2*x2
    endfunction

endmodule

program www_fcode_cn
use GSLD
implicit none
real*16::fn(n),ak(n),x,x1,x2,a,b,a1,b1,a2,b2,answer
integer::i,j,k
call fn0(fn,ak)
a=0
b=1    !积分上下限
a1=0
b1=1
a2=0
b2=1
answer=0
do k=1,n
    do j=1,n
        do i=1,n            answer=answer+ak(i)*ak(j)*ak(k)*y((a+b)/2+(b-a)/2*fn(i),(a1+b1)/2+(b1-a1)/2*fn(j),(a1+b1)/2+(b1-a1)/2*fn(k))
        enddo
    enddo
enddo
answer=answer*(b-a)/2*(b1-a1)/2*(b2-a2)/2
write(*,*)answer
pause
end program www_fcode_cn

pasuka

lz文献查阅力度还不够，对于振荡函数的积分请用自适应高斯-克罗朗德积分一般来说21个点的高斯-克罗朗德积分就够用了，可以先用matlab的quadgk函数体验一下

岸边的鱼

致谢：以下是另外一种求解高斯点和权重系数的程序，由可爱的群友‘红富士’童鞋分享，由我负责整理，标注，悲催的是这个程序里面核心的求解部分我没看懂。希望有能力看懂的群友们，不吝赐教，可以指明里面到底是怎么求解的，不胜感激。

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

!x(n)高斯零点;w(n)权系数;n高斯点个数(可变)   
program main
  implicit none
  integer,parameter:: n=1000   !n高斯点个数(可变)
  integer,parameter::nn=8
  real(kind=nn) ::x1,x2,x(n),w(n)
  real(kind=nn) ::y,f,a
  integer ::i
!------------------------------------------------
!x1,x2积分上下限;变换积分上下限为[-1,1 ]
    x1=-10.0e1
    x2=10.0e1
    call gauleg(-1.d0,1.d0,x,w,n)    !调用求解高斯点程序    
       y=0.0                         !积分结果赋初始值 
    do i=1,n
        y=y+0.5*(x2-x1)*w(i)*f(0.5*(x1+x2)+0.5*(x2-x1)*x(i))                                     !高斯求积公式，已经对被积函数做变量代换处理，
                                     !将被积分区间换到【-1,1】
    end do
    write(*,*)    '积分理论值',sqrt(3.1415926)*exp(-1.d0/4.d0)  !测试计算理论值
    write(*,*)    '数值计算值',y                                !测试数值计算值
    end 
    
function f(x)
    implicit  none
    integer,parameter::nn=8
    real(kind=nn):: f,x
    
       f=exp(-x**2)*cos(x)
       
end function
    
SUBROUTINE gauleg(x1,x2,x,w,n)   !gauleg(-1.0,1.0,x,w,n)
!这x1，x2的传不传的没意思吧，勒让德求0点不就是在[-1,1]上吗？
  implicit none
  integer::n    !计算的阶数
  integer,parameter::nn=8
  real(kind=nn)::x1,x2,x(n),w(n)!积分上下限传递进来
  real(kind=nn),parameter::EPS=3.d-14!系数为何取3，真心诡异
  integer::i,j,m   
  real(kind=8)::p1,p2,p3,pp,xl,xm,z,z1
      m=(n+1)/2          !截取一半，这里取整，这个很厉害
      xm=0.5d0*(x2+x1)     
      xl=0.5d0*(x2-x1)   !将求解区间截取一半，在一半区间上求解高斯点
      do 12 i=1,m        !只计算一半，这个牛逼[0,1]区间求解
         z=cos(3.141592654d0*(i-0.25d0)/(n+0.5d0))   !这是尼玛到底是啥？？
1        continue
          p1=1.0
          p2=0.0
          do 11 j=1,n
            p3=p2
            p2=p1
            p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
11        continue
          pp=n*(z*p1-p2)/(z*z-1.0)   !代替n阶勒让德多项式的导数？？？
          z1=z    
          z=z1-p1/pp
        if(abs(z-z1).gt.EPS)goto 1
        x(i)=xm-xl*z              !代替高斯点
        x(n+1-i)=xm+xl*z          !代替另外一半的高斯点

        w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)  !权重的表达式
        w(n+1-i)=w(i)                    !对称的另一半的权重
    
12    continue
   print*, x,  w    
      return
      
      END

岸边的鱼

pasuka 发表于 2014-6-16 18:08
lz文献查阅力度还不够，对于振荡函数的积分请用自适应高斯-克罗朗德积分一般来说21个点的高斯-克罗朗德积分 ...

终于看懂你的lz是啥意思了，原来是楼主啊，自适应法不知道是不是就是自动变步长法？我其实是有使用的，不要说是21个点，我使用的是100个点的，效果也不是很好，不知道是我用的不好还是什么原因，希望能加个好友再向你请教下。添加方式我已经发你站内信了，谢谢

pasuka

岸边的鱼 发表于 2014-6-17 19:41
终于看懂你的lz是啥意思了，原来是楼主啊，自适应法不知道是不是就是自动变步长法？我其实是有使用的，不 ...

抱歉，单位网络没法上qq
至于gauss-kronrod的基本原理，lz可以自己网上搜索或者查看matlab的quadgk函数代码，会有非常详细的介绍

mangix2010

不明觉厉，收藏了~~

岸边的鱼

pasuka 发表于 2014-6-18 09:41
抱歉，单位网络没法上qq
至于gauss-kronrod的基本原理，lz可以自己网上搜索或者查看matlab的quadgk函数代 ...

好吧，我有时间在学习下matlab吧

岸边的鱼

不会添加图片，我重新编辑

瑶远梦想

个人觉得用迭代改变步长时解决震荡函数积分的好方法！

岸边的鱼

瑶远梦想 发表于 2014-8-3 22:18
个人觉得用迭代改变步长时解决震荡函数积分的好方法！

恩，有道理，现在已经改用自定变步长方法计算了