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

3.3边值问题的直接积分

这节编写时用到了已知结果方程来得到第二个数值，方法不严谨

program main
    implicit none
    integer(8)::n,i,j,k
    real(8),external::su,exact
    real(8)::h,sp(1000),spi(1000),r(1000),haha(1000),m
    open(unit=99,file="xable.csv")
    open(unit=100,file="yable.csv")
    open(unit=101,file="exaction.csv")
    n=500
    h=0.1
    sp(1)=0
    sp(2)=1-0.5*(h+2)*exp(-h)
    sp(2)=0.995*sp(2)
    do j=1,n+1,1
        r(j)=(j-1)*h
        write(99,*)r(j)
    end do
    do i=2,n,1
        sp(i+1)=2*sp(i)-sp(i-1)&        
        +h*h/12*(su(r(i+1))+su(r(i-1))+10*su(r(i)))
    end do
    !!linear correction
    m=(sp(n)-sp(n-100))/(r(n)-r(n-100))
    do k=1,n+1,1
        spi(k)=sp(k)-m*(k-1)*h
        write(100,*)spi(k)
    end do
    do i=1,n+1,1
        haha(i)=exact(r(i))
        write(101,*)haha(i)
    end do
end program main
!!ϕ(r) 满足的 Poisson 方程
function su(x)
    implicit none
    real(8)::x
    real(8)::su
    su=-0.5*x*dexp(-x)
    return 
end 
function exact(y)
    implicit none
    real(8)::y
    real(8)::exact
    exact=1-0.5*(y+2)*dexp(-y)
    return 
end