!########################################################### ! page.53 表3.2 ! ガウス・ルジャンドル積分公式のプログラム ! ! f(x)=1.d0/sqrt(1.d0-x**2) !########################################################### PROGRAM prog08 IMPLICIT REAL(8)(a-h,o-z) IMPLICIT INTEGER(i-n) INTEGER,PARAMETER :: NMAX=200 REAL(8),DIMENSION(NMAX) :: xx2,ww2 REAL(8),DIMENSION(0:NMAX) :: P REAL(8) :: I3 pi=4.d0*datan(1.d0) !******************** DO N=50,200,50 !******************** xx2=0.d0;ww2=0.d0 !**********************************************! ! ガウス・ルジャンドル公式 !**********************************************! CALL ZEROTEN(N,NMAX,xx2) DO i=1,N a=0.d0 CALL legendre(N,xx2(i),P,NMAX) DO k=0,N-1 a=a+dble(2*k+1)*(P(k))**2 END DO ww2(i)=2.d0/a END DO I3=0.d0 DO i=1,N I3=I3+ww2(i)*func(xx2(i)) END DO !******************** E3=dabs(I3-pi) WRITE(*,"('N=',I3,1p3E16.8)") N,E3 !******************** END DO !******************** STOP END PROGRAM prog08 ! ****************************************************************** ! ! ルジャンドル多項式 ! P_j(x) (j=0,1,2,...,N) ! ! ****************************************************************** SUBROUTINE legendre(N,x,P,NMAX) IMPLICIT REAL(8)(a-h,o-z) IMPLICIT INTEGER(i-n) INTEGER,INTENT(IN) :: N,NMAX REAL(8),INTENT(IN) :: x REAL(8),INTENT(OUT),DIMENSION(0:NMAX) :: P P=0.d0 P(0)=1.d0 P(1)=x DO j=1,N-1 P(j+1)=(dble(2*j+1)*x*P(j)-dble(j)*P(j-1))/dble(j+1) END DO RETURN END SUBROUTINE legendre !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo! SUBROUTINE ZEROTEN(N,NMAX,xx) IMPLICIT REAL(8)(a-h,o-z) IMPLICIT INTEGER(i-n) INTEGER,INTENT(IN) :: N,NMAX REAL(8),INTENT(OUT),DIMENSION(NMAX) :: xx REAL(8),DIMENSION(0:NMAX) :: P pi=4.d0*datan(1.d0) xx=0.d0 delta=1.d-10 eps=1.d-10 M=N/2 DO i=1,M a=dsin(pi*dble(N-2*i)/dble(2*N)) b=dsin(pi*dble(N-2*i+2)/dble(2*N)) CALL legendre(N,a,P,NMAX) fa=P(N) CALL legendre(N,b,P,NMAX) fb=P(N) IF(fa*fb>0) THEN WRITE(*,*) 'Error a,b, i==',i ENDIF fc=1.d0 DO WHILE(dabs(fc)>delta.and.dabs(a-b)>eps) diff=dabs(a-b) c=(a+b)/2.d0 CALL legendre(N,c,P,NMAX) fc=P(N) IF(fc*fa>0.d0)THEN a=c ; fa=fc ELSE b=c ; fb=fc END IF END DO 100 CONTINUE xx(N-i+1)=c xx(i)=-c END DO RETURN END SUBROUTINE ZEROTEN !oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FUNCTION func(x) REAL(8),INTENT(IN) :: x REAL(8) :: func func=1.d0/dsqrt(1.d0-x**2) RETURN END FUNCTION func !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!