回 帖 发 新 帖 刷新版面

主题:高斯数值积分如何求被积函数是绝对值的函数吗?

想用高斯数值积分方法求被积函数为绝对值函数的一维积分(调用imsl函数库可以,但是我编的程序用此方法非常不便) 请教各位大虾帮帮忙啊[em1]

回复列表 (共8个回复)

沙发

调用库函数会比你自己编程更麻烦?我十分怀疑。  你若是简单用 Gauss Quadrature 就是一个循环就可以实现,问题是精度你能保证吗?

板凳

因为求这个积分只是我的程序中的一小部分,是作为子程序来调用的,而且其中还涉及到一些公用变量,单纯地用imsl函数库求积分是没问题的, 只不过在整个程序中稍显复杂了些,我只是想简单的用循环来搞定,只不过精度有点问题,可能是积分点的选取问题,但对绝对值函数的高斯积分点怎么选取啊,还请指教下![em18]

3 楼

Gauss Quadrature 的积分节点、权重,只是和权函数、积分区间相关,与被积函数没有什么关系。对某些有奇异点的函数,需要用特殊的算法。 绝对值函数,这个并不能说明它有什么奇异性,比如函数
f(x) = |x|,  这个函数就很普通呀,积分很好处理。

4 楼

例如:一重积分,积分下限是-1.0,上限是1.0,被积函数为f(x) = |x|,积分点和权系数怎么取?
我直接用f(x) =x取积分点和权系数但结果不对,我知道可能是积分点和权系数取的不对,希望能解答下,谢谢![em1]

5 楼

如下积分程序,示范 Legendre-Gauss 积分,但精度不能保证。

! GaussQuadrature
 
 ! written by Zeng Zhuo-Quan
 ! date: 2010 - 06 - 13
 
 Module GaussQuadrature
   implicit none
   
   integer, parameter:: DP = 8
   ! Legendre-Gauss Quadrature [-1, 1]
   integer, parameter:: ND = 15
   real(kind = DP), save:: QX(ND) = (/ &
 -0.9879925180204860D+00, -0.9372733924007061D+00, &
 -0.8482065834104269D+00, -0.7244177313601700D+00, &
 -0.5709721726085385D+00, -0.3941513470775634D+00, &
 -0.2011940939974347D+00,  0.0D0, &
  0.2011940939974347D+00,  0.3941513470775637D+00, &
  0.5709721726085392D+00,  0.7244177313601701D+00, &
  0.8482065834104270D+00,  0.9372733924007057D+00, &
  0.9879925180204853D+00 /)
  
    real(kind = DP), save:: QW(ND) = (/ &
  0.3075324199611685D-01, 0.7036604748810857D-01, &
  0.1071592204671715D+00, 0.1395706779261550D+00, &
  0.1662692058169940D+00, 0.1861610000155613D+00, &
  0.1984314853271117D+00, 0.2025782419255612D+00, &
  0.1984314853271111D+00, 0.1861610000155624D+00, &
  0.1662692058169943D+00, 0.1395706779261532D+00, &
  0.1071592204671729D+00, 0.7036604748810787D-01, &
  0.3075324199611745D-01 /)
  
contains 
  subroutine Gass_Quad(pf, a, b, ans)
    implicit none
!    real(kind = DP), external:: pf   ! procedure(f):: pf
    procedure(f), pointer:: pf   ! Fortran 2003 standard
    real(kind = DP), intent(in):: a   ! lower bound
    real(kind = DP), intent(in):: b   ! upper bound
    real(kind = DP), intent(out):: ans
    !
    integer:: i   ! for loop
    real(kind = DP):: apb, amb
    real(kind = DP):: xi, wi
    
    apb = 0.5D0 * (a + b)
    amb = 0.5D0 * (a - b)
    ans = 0.0D0
    Do i = 1, ND, 1
      xi = apb - amb * QX(i)
      wi = -amb * QW(i)
      ans = ans + wi * pf(xi)
    End Do   ! i
    return
  end subroutine
  
  real(kind = DP) function f(x)
    implicit none
    real(kind = DP), intent(in):: x
    f = abs(x)
    return
  end function 
  
End Module


program main 
  use GaussQuadrature
  implicit none
  real(kind = DP):: a, b
  real(kind = DP):: ans  
  procedure(f), pointer:: pf
  a = -1.0D0   ! lower bound
  b = 1.0D0    ! upper bound
  
  pf => f   ! Fortran 2003 standard
  call Gass_Quad(pf, a, b, ans)
! call Gass_Quad(f, a, b, ans)
  write(*, *) "Quadrature result: ", ans
  
  stop
end program main

6 楼

想给你继续评分来着 可系统不让整不好意思了!这里感谢先!!希望以后有问题继续向你请教。

7 楼

你好,我也是在用fortran编程时用到了三维积分的问题,自己只用简单的三重循环实现高斯积分最后不同点结果相差很大,是不是自己编程积分容易发散,是不是应该用徐士良的算法集里面的多重积分代码,可是看不大懂,能不能解释一下源代码啊?谢谢,非常非常感谢
SUBROUTINE FGAUS(N,JS,X,FS,F,S)
    DIMENSION JS(N),X(N)
    DIMENSION T(5),C(5),D(2,11),CC(11),IS(2,11)
    DOUBLE PRECISION X,F,S,T,C,D,CC,DN,UP,P
    DATA T/-0.9061798459,-0.5384693101,0.0,
     *         0.5384693101,0.9061798459/
    DATA C/0.2369268851,0.4786286705,0.5688888889,
     *         0.4786286705,0.2369268851/
    M=1
    D(1,N+1)=1.0
    D(2,N+1)=1.0
10    DO 20 J=M,N
      CALL FS(J,N,X,DN,UP)
      D(1,J)=0.5*(UP-DN)/JS(J)
      CC(J)=D(1,J)+DN
      X(J)=D(1,J)*T(1)+CC(J)
      D(2,J)=0.0
      IS(1,J)=1
      IS(2,J)=1
20    CONTINUE

    J=N
30    K=IS(1,J)
    IF (J.EQ.N) THEN
      P=F(N,X)
    ELSE
      P=1.0
    END IF
    D(2,J)=D(2,J+1)*D(1,J+1)*P*C(K)+D(2,J)
    IS(1,J)=IS(1,J)+1
    IF (IS(1,J).GT.5) THEN
      IF (IS(2,J).GE.JS(J)) THEN
        J=J-1
        IF (J.EQ.0) THEN
          S=D(2,1)*D(1,1)
          RETURN
        END IF
        GOTO 30
      END IF
      IS(2,J)=IS(2,J)+1
      CC(J)=CC(J)+D(1,J)*2.0
      IS(1,J)=1
    END IF
    K=IS(1,J)
    X(J)=D(1,J)*T(K)+CC(J)
    IF (J.EQ.N) GOTO 30
    M=J+1
    GOTO 10
    END

8 楼

把这个绝对值函数化成分段函数啊,-1<x<0,f=-x;1>x>0,f=x

我来回复

您尚未登录,请登录后再回复。点此登录或注册