主题:高斯数值积分如何求被积函数是绝对值的函数吗?
kings
[专家分:0] 发布于 2010-06-05 16:42:00
想用高斯数值积分方法求被积函数为绝对值函数的一维积分(调用imsl函数库可以,但是我编的程序用此方法非常不便) 请教各位大虾帮帮忙啊[em1]
回复列表 (共8个回复)
沙发
asymptotic [专家分:16630] 发布于 2010-06-09 17:06:00
调用库函数会比你自己编程更麻烦?我十分怀疑。 你若是简单用 Gauss Quadrature 就是一个循环就可以实现,问题是精度你能保证吗?
板凳
kings [专家分:0] 发布于 2010-06-10 10:23:00
因为求这个积分只是我的程序中的一小部分,是作为子程序来调用的,而且其中还涉及到一些公用变量,单纯地用imsl函数库求积分是没问题的, 只不过在整个程序中稍显复杂了些,我只是想简单的用循环来搞定,只不过精度有点问题,可能是积分点的选取问题,但对绝对值函数的高斯积分点怎么选取啊,还请指教下![em18]
3 楼
asymptotic [专家分:16630] 发布于 2010-06-10 10:47:00
Gauss Quadrature 的积分节点、权重,只是和权函数、积分区间相关,与被积函数没有什么关系。对某些有奇异点的函数,需要用特殊的算法。 绝对值函数,这个并不能说明它有什么奇异性,比如函数
f(x) = |x|, 这个函数就很普通呀,积分很好处理。
4 楼
kings [专家分:0] 发布于 2010-06-13 15:42:00
例如:一重积分,积分下限是-1.0,上限是1.0,被积函数为f(x) = |x|,积分点和权系数怎么取?
我直接用f(x) =x取积分点和权系数但结果不对,我知道可能是积分点和权系数取的不对,希望能解答下,谢谢![em1]
5 楼
asymptotic [专家分:16630] 发布于 2010-06-13 17:56:00
如下积分程序,示范 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 楼
kings [专家分:0] 发布于 2010-06-13 18:56:00
想给你继续评分来着 可系统不让整不好意思了!这里感谢先!!希望以后有问题继续向你请教。
7 楼
pingnuaa [专家分:0] 发布于 2011-12-08 11:08:00
你好,我也是在用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 楼
alsoran [专家分:760] 发布于 2011-12-08 11:40:00
把这个绝对值函数化成分段函数啊,-1<x<0,f=-x;1>x>0,f=x
我来回复