FORTRAN+计算物理学学习日记（12）

第四章 特殊函数和高斯求积法

legendre多项式

program main
    implicit none
    integer(8)::p,n,i,j,m
    real(8)::h,x(3000),q(3000,3000),tol,xx(300)
    real(8),external::legendre
    p=4
    n=1000
    h=2.0/n
    tol=1e-10
    m=10
    open(unit=100,file="xable.csv")
    open(unit=101,file="y1.csv")
    open(unit=102,file="y2.csv")
    open(unit=103,file="y3.csv")
    open(unit=104,file="y4.csv")
    open(unit=105,file="y5.csv")
    open(unit=106,file="zeros.csv")
    do j=1,p+1,1
        do i=1,n+1,1
            x(i)=-1.0+(i-1)*h
            q(i,j)=legendre(j-1,x(i))
        end do
    end do
    do i=1,n+1,1
        write(100,*)x(i)
        write(101,*) q(i,1)
        write(102,*) q(i,2)
        write(103,*) q(i,3)
        write(104,*) q(i,4)
        write(105,*) q(i,5)
    end do
    close(100)
    close(101)
    close(102)
    close(103)
    close(104)
    close(105)
    call zeros(m,tol,xx)
    do i=1,m,1
        write(106,*)xx(i)
    end do
    close(106)
end program main
function legendre(lamx,x)
    real(8)::x,p(300)
    integer(8)::lamx,i
    real(8)::legendre
    p(1)=1.0
    p(2)=x
    do i=2,lamx,1
        p(i+1)=((2.0*i-1.0)*x*p(i)-(i-1)*p(i-1))/i
    end do
    legendre=p(lamx+1)
    return 
end 
subroutine zeros(m,tol,x)
  implicit none
  real(8),intent(in)::tol
  integer(8),intent(in)::m
  real(8)::a,dx,c,b
  real(8),intent(out)::x(300)
  integer(8)::i
  real(8),external::legendre
  a=-1.0
  do i=1,m,1
      dx=0.01
      a=a+dx
      b=legendre(m,a)
      do while(abs(dx)>tol)
          a=a+dx
          c=legendre(m,a)
          if(b*c>0)then
          else
              a=a-dx
              dx=dx/2.0
          end if
      end do
      x(i)=a
  end do
  return
end subroutine zeros

 laguerre多项式求积法

program main
    implicit none
    integer(8)::p,n,i,j,m
    real(8)::h,x(3000),q(3000,3000),tol,xx(3000)
    real(8),external::laguerre
    p=4
    n=100
    h=5.0/n
    tol=1e-10
    m=10
    open(unit=100,file="xable.csv")
    open(unit=101,file="y1.csv")
    open(unit=102,file="y2.csv")
    open(unit=103,file="y3.csv")
    open(unit=104,file="y4.csv")
    open(unit=105,file="y5.csv")
    open(unit=106,file="zeros.csv")
    do j=1,p+1,1
        do i=1,n+1,1
            x(i)=(i-1)*h
            q(i,j)=laguerre(j-1,x(i))
        end do
    end do
    do i=1,n+1,1
        write(100,*) x(i)
        write(101,*) q(i,1)
        write(102,*) q(i,2)
        write(103,*) q(i,3)
        write(104,*) q(i,4)
        write(105,*) q(i,5)
    end do
    close(100)
    close(101)
    close(102)
    close(103)
    close(104)
    close(105)
    call zeros(m,tol,xx)
    do i=1,m,1
        write(106,*)xx(i)
    end do
    close(106)
end program main
function laguerre(nmax,x)
  implicit none
  real(8)::x,L(3000)
  integer(8)::nmax,i
  real(8)::laguerre
  L(1)=1.0
  L(2)=1.0-x
  do i=2,nmax,1
      L(i+1)=((2.0*i-1.0-x)*L(i)-(i-1.0)*L(i-1))/i
  end do
  laguerre=L(nmax+1)
  return 
end 
subroutine zeros(n,tol,x)
  implicit none
  integer(8),intent(in)::n
  integer(8)::i
  real(8),intent(in)::tol
  real(8),intent(out)::x(3000)
  real(8)::a,b,dx,c
  real(8),external::laguerre
  a=0.0
  do i=1,n,1
      dx=0.01
      a=a+dx
      b=laguerre(n,a)
      do while(abs(dx)>tol)
          a=a+dx
          c=laguerre(n,a)
          if(b*c>0)then
          else
              a=a-dx
              dx=dx/2.0
          end if
      end do
      x(i)=a
  end do
  return
end subroutine zeros

hermite求积法

program main
    implicit none
    integer(8)::n,m,i,j,z
    real(8)::b,h,x(1000),q(1000,1000),tol,xx(1000)
    real(8),external::hermite
    n=5
    m=200
    b=2.0
    z=10
    h=2.0*b/m
    tol=1e-10
    open(unit=100,file="xable.csv")
    open(unit=101,file="y1.csv")
    open(unit=102,file="y2.csv")
    open(unit=103,file="y3.csv")
    open(unit=104,file="y4.csv")
    open(unit=105,file="y5.csv")
    open(unit=106,file="y6.csv")
    open(unit=107,file="zeros.csv")
    do j=1,n+1,1
        do i=1,m+1,1
            x(i)=-b+(i-1)*h
            q(i,j)=hermite(j-1,x(i))
        end do
    end do
    call zeros(z,tol,xx)
    do i=1,z,1
        write(107,*)xx(i)
    end do
    do i=1,m+1,1
        write(100,*)x(i)
        write(101,*)q(i,1)
        write(102,*)q(i,2)
        write(103,*)q(i,3)
        write(104,*)q(i,4)
        write(105,*)q(i,5)
        write(106,*)q(i,6)
    end do
    close(100)
    close(101)
    close(102)
    close(103)
    close(104)
    close(105)
    close(106)
    close(107)
end program main
function hermite(nmax,x)
    real(8)::x,h(1000)
    integer(8)::nmax,i
    real(8)::hermite
    h(1)=1.0
    h(2)=2.0*x
    do i=2,nmax,1
        h(i+1)=2.0*x*h(i)-2.0*(i-1)*h(i-1)
    end do
    hermite=h(nmax+1)
end 
subroutine zeros(n,tol,xx)
    implicit none
    integer(8),intent(in)::n
    real(8),intent(in)::tol
    real(8),intent(out)::xx(1000)
    real(8)::a,dx,c,b
    real(8),external::hermite
    integer(8)::i
    a=-n+0.0
    do i=1,n,1
        dx=0.01
        a=a+dx
        b=hermite(n,a)
        do while(abs(dx)>tol)
            a=a+dx
            c=hermite(n,a)
            if(b*c>0)then
            else
                a=a-dx
                dx=dx/2.0
            end if
        end do
        xx(i)=a
    end do
end subroutine zeros

 

Gauss-Hermite 求积公式

program main
    implicit none
    integer(8)::n,k,j
    real(8)::z,x(1000),tol,int,m,pi,A(1000)
    real(8),external::f3,hermite
    n=10
    tol=1e-15
    m=1.0
    int=0.0
    pi=3.1415926
    open(unit=100,file="result.csv")
    call zeros(n,tol,x)
    do j=1,n,1
        m=m*j
    end do
    do k=1,n,1
       A(k)=2**(n-1)*m*dsqrt(pi)/((n*hermite(n-1,x(k)))*(n*hermite(n-1,x(k))))
        int=int+A(k)*f3(x(k))
    end do
    write(100,*)int
    close(100)
end program main
function f3(x)
    implicit none
    real(8)::x
    real(8)::f3
    f3=cos(x)
end 
function hermite(nmax,x)
    real(8)::x,h(1000)
    integer(8)::nmax,i
    real(8)::hermite
    h(1)=1.0
    h(2)=2.0*x
    do i=2,nmax,1
        h(i+1)=2.0*x*h(i)-2.0*(i-1)*h(i-1)
    end do
    hermite=h(nmax+1)
end 
subroutine zeros(n,tol,xx)
    implicit none
    integer(8),intent(in)::n
    real(8),intent(in)::tol
    real(8),intent(out)::xx(1000)
    real(8)::a,dx,c,b
    real(8),external::hermite
    integer(8)::i
    a=-n+0.0
    do i=1,n,1
        dx=0.01
        a=a+dx
        b=hermite(n,a)
        do while(abs(dx)>tol)
            a=a+dx
            c=hermite(n,a)
            if(b*c>0)then
            else
                a=a-dx
                dx=dx/2.0
            end if
        end do
        xx(i)=a
    end do
end subroutine zeros    

Gauss-Laguerre 求积公式

program main
    implicit none
    integer(8)::n,k
    real(8)::x(3000),int,A(3000),tol
    real(8),external::laguerre,f2
    n=5
    tol=1e-10
    call zeros(n,tol,x)
    int=0.0
    open(unit=100,file="result.csv")
    do k=1,n,1
        A(k)=x(k)/((n*laguerre(n-1,x(k)))*(n*laguerre(n-1,x(k))))
        int=int+A(k)*f2(x(k))
    end do
    write(100,*)int
    close(100)
end program main
function f2(x)
    implicit none
    real(8)::x
    real(8)::f2
    f2=x*x*x+x*x
    return
end 
function laguerre(nmax,x)
  implicit none
  real(8)::x,L(3000)
  integer(8)::nmax,i
  real(8)::laguerre
  L(1)=1.0
  L(2)=1.0-x
  do i=2,nmax,1
      L(i+1)=((2.0*i-1.0-x)*L(i)-(i-1.0)*L(i-1))/i
  end do
  laguerre=L(nmax+1)
  return 
end 
subroutine zeros(n,tol,x)
  implicit none
  integer(8),intent(in)::n
  integer(8)::i
  real(8),intent(in)::tol
  real(8),intent(out)::x(3000)
  real(8)::a,b,dx,c
  real(8),external::laguerre
  a=0.0
  do i=1,n,1
      dx=0.01
      a=a+dx
      b=laguerre(n,a)
      do while(abs(dx)>tol)
          a=a+dx
          c=laguerre(n,a)
          if(b*c>0)then
          else
              a=a-dx
              dx=dx/2.0
          end if
      end do
      x(i)=a
  end do
  return
end subroutine zeros    

Gauss-Legendre 求积公式

program main
    implicit none
    integer(8)::n,k
    real(8)::x(1000),int,A(1000),h
    real(8),external::legendre,f1
    n=10
    h=1e-15
    call zeros(n,h,x)
    int=0.0
    open(unit=100,file="result.csv")
    do k=1,n,1
        A(k)=2.0*(1.0-x(k)*x(k))/((n*legendre(n-1,x(k)))*(n*legendre(n-1,x(k))))
        int=int+A(k)*f1(x(k))
    end do
    write(100,*)int
end program main
function f1(x)
    implicit none
    real(8)::x
    real(8)::f1
    f1=3.0/2.0*sqrt(3.0/2.0*x+5.0/2.0)
end 
function legendre(lamx,x)
    real(8)::x,p(300)
    integer(8)::lamx,i
    real(8)::legendre
    p(1)=1.0
    p(2)=x
    do i=2,lamx,1
        p(i+1)=((2.0*i-1.0)*x*p(i)-(i-1)*p(i-1))/i
    end do
    legendre=p(lamx+1)
    return 
end 
subroutine zeros(m,tol,x)
  implicit none
  real(8),intent(in)::tol
  integer(8),intent(in)::m
  real(8)::a,dx,c,b
  real(8),intent(out)::x(300)
  integer(8)::i
  real(8),external::legendre
  a=-1.0
  do i=1,m,1
      dx=0.01
      a=a+dx
      b=legendre(m,a)
      do while(abs(dx)>tol)
          a=a+dx
          c=legendre(m,a)
          if(b*c>0)then
          else
              a=a-dx
              dx=dx/2.0
          end if
      end do
      x(i)=a
  end do
  return
end subroutine zeros