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

3.5打靶法

边值问题

program main
    implicit none
    integer(8)::n,i
    real(8)::h,error,delta,ddelta,xs(100),yr(100),yr1(100),ka(100),kb(100),yn(100),ynn(100),re(100),pi
    real(8),external:: g5
    pi=3.1415926
    n=40
    h=1.0/n
    error=10e-5
    delta=0.8
    ddelta=0.1
    open(unit=100,file="xable.csv")
    open(unit=101,file="numerov.csv")
    open(unit=102,file="exact.csv")
    do i=1,n+1,1
        xs(i)=h*(i-1)
        write(100,*)xs(i)
    end do
    close(100)
    !! The Numerov Algorithm
    call interp(delta,h,n,ynn)
    do while(ddelta>error)
        delta=delta-ddelta
        call interp(delta,h,n,yn)
        if((yn(n+1)-1)*(ynn(n+1)-1)>0)then
        else
            delta=delta+ddelta
            ddelta=ddelta/2.0
        end if
    end do
    do i=1,n+1,1
        write(101,*)yn(i)
    end do
    close(101)
    !!The exact result
    do i=1,n+1,1
        re(i)=g5(xs(i))
        write(102,*)re(i)
    end do
    close(102)
    end program
    
function g5(x)
    implicit none
    real(8)::x,pi
    real(8)::g5
    pi=3.1415926
    g5=cos(1.0/2*pi*x)+2.0*sin(1.0/2*pi*x)-1.0
    return
end    
subroutine interp(delta,h,n,yx)
    integer(8)::i
    real(8)::con,pi
    real(8),intent(in)::delta
    real(8),intent(in)::h
    integer(8),intent(in)::n
    real(8),intent(out)::yx(100)
    yx(1)=0
    yx(2)=delta
    pi=3.1415926
    con=h*h*pi*pi/48.0
    do i=2,n,1
        yx(i+1)=2.0*(1.0-5*con)*yx(i)-(1+con)*yx(i-1)-h*h*pi*pi/4.0
        yx(i+1)=yx(i+1)/(1+con)
    end do
    return
end subroutine interp

本征值问题

program main
    implicit none
    integer(8)::i,n
    real(8)::h,delta,tol,k,dk,x(500),y(500),p1(500),p2(500),a,pp(500),p(100)
    real(8),external::g4
    n=50
    h=1.0/n
    delta=0.6
    tol=10**(-5)
    k=12.0
    dk=1.0
    open(unit=100,file="xable.csv")
    open(unit=101,file="normalization.csv")
    open(unit=102,file="exact.csv")
    call interp(k,h,n,delta,p1)
    do while(abs(dk)>tol)
        k=k+dk
        call interp(k,h,n,delta,p2)
        if(p1(n+1)*p2(n+1)>0)then
        else
            k=k-dk
            dk=dk/2.0
        end if
    end do
    do i=1,n+1,1
        x(i)=(i-1)*h
        write(100,*)x(i)
    end do
    !!the normalization
    a=0
    do i=1,n,1
        a=a+h/2*(p2(i)*p2(i)+p2(i+1)*p2(i+1))
    end do
    do i=1,n+1,1
        p2(i)=p2(i)/dsqrt(a)
        write(101,*)p2(i)
    end do
    !!exact solution
    do i=1,n+1,1
        y(i)=g4(x(i))
        write(102,*)y(i)
    end do
    close(100)
    close(101)
    close(102)
    
end program main
    
function g4(x)
    implicit none
    real(8)::x,pi
    real(8)::g4
    pi=3.1415926
    g4=sqrt(2.0)*dsin(4.0*pi*x)
    return
    end
subroutine interp(k,h,n,delta,p)
    implicit none
    real(8),intent(in)::k,h,delta
    integer(8),intent(in)::n
    real(8),intent(out)::p(100)
    real(8) con
    integer(8) j
    p(1)=0
    p(2)=delta
    con=(k*h)*(k*h)/12.0
    do j=2,n,1
        p(j+1)=2.0*(1-5*con)*p(j)-(1+con)*p(j-1)
        p(j+1)=p(j+1)/(1+con)
    end do
end subroutine interp