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

阿玛兰妲.com

已于 2024-08-15 13:17:49 修改

阅读量207

点赞数

文章标签：学习几何学

于 2022-06-11 11:39:00 首次发布

本文链接：https://blog.csdn.net/m0_68623139/article/details/125229625

版权

本文详细介绍了使用Adams两步法、四步法及Adams-Moulton两步法求解常微分方程的程序实现，并通过与精确解的比较，展示了不同算法的精度。重点讨论了Adams算法在数值解微分方程中的应用和步骤。

摘要由CSDN通过智能技术生成

2.3多步法求解常微分方程

program main
    implicit none
    integer(8)::n,m,i,k
    real(8),external::fa,fb
    real(8)::h,x(200),ya(200),yb(200),yc(200,200),yd(200)
    n=100
    m=100
    h=4.0/n
    open(unit=99,file="xable.csv")
    open(unit=100,file="Adams-two.csv")
    open(unit=101,file="Adams-four.csv")
    open(unit=102,file="Adams-Moulton two.csv")
    open(unit=103,file="exact.csv")
    do i=1,n+1,1
        x(i)=(i-1)*h
        write(99,*)x(i)
    end do
    close(99)
    !The Adams two-step algorithm
    ya(1)=1
    ya(2)=1-h*h/2.0
    write(100,*)ya(1)
    write(100,*)ya(2)
    do i=2,n,1
        ya(i+1)=ya(i)+h*(3.0/2.0*fa(x(i-1),ya(i))-fa(x(i-1),ya(i-1))/2.0)
        write(100,*)ya(i+1)
    end do
    close(100)
    !The Adams four-step algorithm
    yb(1)=1
    write(101,*)yb(1)
    do i=1,3,1
        yb(i+1)=yb(i)+h*fa(x(i),yb(i))
        write(101,*)yb(i+1)
    end do
    do i=4,n,1
        yb(i+1)=yb(i)+h/24.0*(55.0*fa(x(i),yb(i))-59.0*fa(x(i-1),yb(i-1))&
        +37.0*fa(x(i-2),yb(i-2))-9.0*fa(x(i-3),yb(i-3)))
        write(101,*)yb(i+1)
    end do
    close(101)
    !The Adams-Moulton two-step algorithm
    yc(1,1)=1
    yc(1,2)=1-h*h/2.0
    do i=2,n,1
        yc(1,i+1)=yc(1,i)+h*(3.0/2.0*fa(x(i),yc(1,i))&
        -fa(x(i-1),yc(1,i-1))/2.0)
    end do
    do i=2,n,1
        do k=1,m,1
            yc(k+1,1)=1
            yc(k+1,2)=yc(1,2)
            yc(k+1,i+1)=yc(1,i)+h/12.0*(5.0*fa(x(i+1),yc(k,i+1))&
            +8.0*fa(x(i),yc(k,i))-9.0*fa(x(i-1),yc(k,i-1)))   
        end do
    end do
    do i=1,n+1,1
        write(102,*)yc(m+1,i)
    end do
    close(102)
    !The exact solution
    do i=1,n+1,1
        yd(i)=fb(x(i))
        write(103,*)yd(i)
    end do
    close(103)
end program main

    !!exact solution 的function函数
function fb(x)
    implicit none
    real(8)::x
    real(8)::fb
    fb=dexp(-x*x/2.0)
    return 
end
!!-x*y的function函数
function fa(x,y)
    implicit none
    real(8)::x,y
    real(8)::fa
    fa=-x*y
    return 
end