我繼續查解答，從...
Numerical Recipes in Fortran 77 Second Edition
4.4 Improper Integrals
裡面提供SUBRUTINE，
我還不知道要怎麼用？
比如說 aa=1, bb=無限大 積分函數是 Exp(-4x)
在主程式要怎麼CALL？
如下
SUBROUTINE midexp(funk,aa,bb,s,n)
! aa 是積分下限，bb 是積分上限
! n 是整數
! s 是實數
! funk 是自訂的函式
INTEGER n
REAL aa,bb,s,funk
EXTERNAL funk
INTEGER it,j
REAL ddel,del,sum,tnm,x,func,a,b
func(x)=funk(-log(x))/x
b=exp(-aa)
a=0.
if (n.eq.1) then
s=(b-a)*func(0.5*(a+b))
else
it=3**(n-2)
tnm=it
del=(b-a)/(3.*tnm)
ddel=del+del
x=a+0.5*del
sum=0.
do j=1,it
sum=sum+func(x)
x=x+ddel
sum=sum+func(x)
x=x+del
enddo
s=(s+(b-a)*sum/tnm)/3.0
endif
RETURN
END
※ 引述《phystw (phys.tw)》之銘言：
: B：積分上限是無限大
: A：積分下限為有限整數假設是0好了
: 請問大家，積分上下限為無限的的狀況該怎麼處理？
: 我想~不應該隨便給一個很大的數值，結果會不一樣。
: 以下我引用彭國倫FORTRAN90的範例
: 積分上限給一個很大的值，
: A=0.0 ! 積分下限
: B=1.0E+06 ! 積分上限
: 積分函數 常數*Exp(-4*x)
: PROGRAM main
: !implicit none
: REAL Pi
: PARAMETER(Pi=3.1415926)
: REAL F, Cross_section_const
: EXTERNAL F, Cross_section_const !補充宣告說明F, Cross_section_const 是函式
: EXTERNAL SIMPSON_INT !補充宣告說明 SIMPSON_INT 是函式
: EXTEREAL A, B ! 積分上限
: REAL ANS ! 積分結果
: A=0.0 ! 積分下限
: B=1.0E+06 ! 積分上限
: ANS = SIMPSON_INT(A, B, F) * &
: SIMPSON_INT(A, B, Cross_section_const)! 常數做積分
: WRITE(*,*) '積分結果:', ANS
: PAUSE
: STOP
: END Program main
: !C
: !C 積分函數
: REAL FUNCTION F(X) ! 自訂函式宣告
: implicit none
: REAL X
: F = Exp(-4*X)
: RETURN
: END
: REAL FUNCTION Cross_section_const ! 自訂函式宣告! 常數做積分
: implicit none
: Cross_section_const = 5.0
: RETURN
: END
: !C
: !C 辛普森法積分函數
: !C
: REAL FUNCTION SIMPSON_INT(A, B, FUNC)
: Implicit None
: REAL A, B
: REAL FUNC
: EXTERNAL FUNC
: INTEGER INTERVALS
: PARAMETER(INTERVALS=10000)
: REAL C
: REAL SUM
: REAL STEP
: REAL STEP2
: STEP = (B-A)/INTERVALS
: STEP2 = STEP*2
: SUM = FUNC(A) + FUNC(B) ! 給一個函數初始值
: DO C = A + STEP, B - STEP, STEP2
: SUM = SUM + 4*FUNC(C)
: EndDo
: DO C = A + STEP2, B - STEP2, STEP2
: SUM = SUM + 2*FUNC(C)
: EndDo
: SIMPSON_INT = SUM*STEP/3.0
: RETURN
: END