[Fortran] 纯文本查看 复制代码
program mytp12
implicit none
integer:: n
double precision,allocatable::xk(:),wk(:)
write(6,*)'input order of quadrature n'
read(5,*)n
allocate(xk(n),wk(n))
call GL_xw(n,xk,wk)
end program mytp12
double precision function leg(n,x,dp)
implicit none
integer :: k,n
double precision:: dp,p(0:n),x
p(0) = 1.d0
p(1) = x
do k = 2, n
p(k)=((2.d0*k-1.d0)*x*p(k-1)-(k-1.d0)*p(k-2))/k
end do
leg=p(n)
dp=(n*p(n-1)-n*x*p(n))/(1-x**2)
end function leg
subroutine GL_xw(n,xk,wk)
implicit none
integer n,k
double precision xk(n),wk(n)
double precision, parameter :: eps = 1.d-15
double precision dp,dx,leg,p,pi,x
pi = acos(-1.d0)
do k = 1, n
x=cos((4.d0*k-1.d0)*pi/(4.d0*n+2.d0))
do
p = leg(x,n,dp)
dx = p / dp
x = x - dx
if (abs (dx) < eps) exit
end do
xk(k) = x
wk(k)=2.d0/((1.d0-x**2)*dp*dp)
end do
do k = 1, n
write (6,'("xk(",i3,")= ",d15.8," wk(",i3,")= ",d15.8)') k,xk(k),k,wk(k)
end do
return
endsubroutine
! --------------------------------------------------
! Silverfrost FTN95 for Microsoft Visual Studio
! Free Format FTN95 Source File
! --------------------------------------------------
