Re: [問題] 關於積分上下限為無限大

看板Fortran作者phystw (phys.tw)時間10年前 (2014/03/14 00:22)推噓0(0推 0噓 0→)

留言0則, 0人參與討論串2/2 (看更多)

我繼續查解答，從... 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 -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 140.113.127.84