5. Ordinary Differential Equation

求解问题

Write a program to solve the following ordinary differential equation by

  • basic Euler method
  • improved Euler method
  • four-order Runge-Kutta method

and calculate y(1.5) with stepsize=0.1, 0.1/2, 0.1/4, 0.1/8

Compare it with analytic solution (in figure)

主程序

program main
    use ODE
    implicit none
    real :: x0=0.0,y0=3.0,a=0.0,b=1.5,h=0.1
    integer :: n,i,j
    character(len=512) :: filename


    print *,'The analytic solution of y(1.5)= ',analytic_solution(1.5)
    print *

    print *,'Basic Euler Method'
    print '(5x,a,f3.1,3x,a,f3.1,3x,a,f3.1,3x,a,f3.1)','Set initial  x0= ',x0,'y0= ',y0,'a= ',a,'b= ',b
    do j=0,3
        h=0.1/2**j
        print '(5x,a,f6.4)','Set h= ',h
        call init(x0,y0,a,b,h)
        n=size(x)-1
        call Basic_Euler
        print '(8x,a,f10.8)','y(1.5)= ',y(n)
        write(filename, *) j
        filename='BEM_'//Trim(AdjustL(filename))//'.txt'
        open(101,file=filename)
        write(101,'(f3.1,f10.8)') (x(i),y(i),i=0,n)
        close(101)
        deallocate(x,y)
    end do
    print *

    print *,'Improved Euler Method'
    print '(5x,a,f3.1,3x,a,f3.1,3x,a,f3.1,3x,a,f3.1)','Set initial  x0= ',x0,'y0= ',y0,'a= ',a,'b= ',b
    do j=0,3
        h=0.1/2**j
        print '(5x,a,f6.4)','Set h= ',h
        call init(x0,y0,a,b,h)
        n=size(x)-1
        call Improved_Euler
        print '(8x,a,f10.8)','y(1.5)= ',y(n)
        write(filename, *) j
        filename='IEM_'//Trim(AdjustL(filename))//'.txt'
        open(101,file=filename)
        write(101,'(f3.1,f10.8)') (x(i),y(i),i=0,n)
        close(101)
        deallocate(x,y)
    end do
    print *

    print *,'Four Order Runge-Kutta Method'
    print '(5x,a,f3.1,3x,a,f3.1,3x,a,f3.1,3x,a,f3.1)','Set initial  x0= ',x0,'y0= ',y0,'a= ',a,'b= ',b
    do j=0,3
        h=0.1/2**j
        print '(5x,a,f6.4)','Set h= ',h
        call init(x0,y0,a,b,h)
        n=size(x)-1
        call Four_Order_Runge_Kutta
        print '(8x,a,f10.8)','y(1.5)= ',y(n)
        write(filename, *) j
        filename='RKM_'//Trim(AdjustL(filename))//'.txt'
        open(101,file=filename)
        write(101,'(f3.1,f10.8)') (x(i),y(i),i=0,n)
        close(101)
        deallocate(x,y)
    end do
    print *

end program main

求解ODE方程的关键方法写在ODE模块中：

module ODE
    implicit none
    private
    public :: x,y,init,Basic_Euler,Improved_Euler,Four_Order_Runge_Kutta,analytic_solution

    real,allocatable :: x(:),y(:)
    real :: x0,y0,a,b,h
    integer :: n,i

contains

    function f(x,y)
        implicit none
        real :: f,x,y
        f=-x*x*y*y
    end function f

    function analytic_solution(x) result(f)
        implicit none
        real :: f,x
        f=3.0/(1+x*x*x)
    end function

    subroutine init(x0_,y0_,a_,b_,h_)
        implicit none
        real :: x0_,y0_,a_,b_,h_
        x0=x0_
        y0=y0_
        a=a_
        b=b_
        h=h_
        n=int((b-a)/h)
        allocate(x(0:n),y(0:n))
        x=(/ (a+i*h,i=0,n) /)
        y=0
        y(0)=y0
    end subroutine init

    subroutine Basic_Euler()
        implicit none
        do i=1,n
            y(i)=y(i-1)+h*f(x(i-1),y(i-1))
        end do
    end subroutine Basic_Euler

    subroutine Improved_Euler()
        implicit none
        real :: y_
        do i=1,n
            y_=y(i-1)+h*f(x(i-1),y(i-1))
            y(i)=y(i-1)+h/2*(f(x(i-1),y(i-1))+f(x(i),y_))
        end do
    end subroutine Improved_Euler

    subroutine Four_Order_Runge_Kutta()
        implicit none
        real :: k1,k2,k3,k4
        do i=1,n
            k1=f(x(i-1),y(i-1))
            k2=f(x(i-1)+h/2,y(i-1)+h/2*k1)
            k3=f(x(i-1)+h/2,y(i-1)+h/2*k2)
            k4=f(x(i-1)+h,y(i-1)+h*k3)
            y(i)=y(i-1)+h/6*(k1+2*k2+2*k3+k4)
        end do
    end subroutine Four_Order_Runge_Kutta
end module ODE

其中init方法是用来初始化：

    subroutine init(x0_,y0_,a_,b_,h_)
        implicit none
        real :: x0_,y0_,a_,b_,h_
        x0=x0_
        y0=y0_
        a=a_
        b=b_
        h=h_
        n=int((b-a)/h)
        allocate(x(0:n),y(0:n))
        x=(/ (a+i*h,i=0,n) /)
        y=0
        y(0)=y0
    end subroutine init

analytic_solution是数值解：

    function analytic_solution(x) result(f)
        implicit none
        real :: f,x
        f=3.0/(1+x*x*x)
    end function

Basic Euler Method

流程图：


原理：
![][3]
[3]: http://latex.codecogs.com/gif.latex?\begin{cases}%20y(x_{n+1})=y(x_n)+hf(x_n,y(x_n))+O(h^2)%20\%20y(x_0)=y_0%20\end{cases}

代码：

    subroutine Basic_Euler()
        implicit none
        do i=1,n
            y(i)=y(i-1)+h*f(x(i-1),y(i-1))
        end do
    end subroutine Basic_Euler

输出结果：


Basic Euler 方法

图例：


不同h每步迭代结果

Improved Euler Method

流程图：


原理：
![][4]
[4]: http://latex.codecogs.com/gif.latex?\begin{cases}%20\bar%20y_{n+1}=y_n+hf(x_n,y_n)%20\%20y(x_{n+1})=y(x_n)+\dfrac{h}{2}[f(x_n,y(x_n))+f(x_{n+1},\bar%20y_{n+1})]+O(h^2)%20\%20y(x_0)=y_0%20\end{cases}

代码：

    subroutine Improved_Euler()
        implicit none
        real :: y_
        do i=1,n
            y_=y(i-1)+h*f(x(i-1),y(i-1))
            y(i)=y(i-1)+h/2*(f(x(i-1),y(i-1))+f(x(i),y_))
        end do
    end subroutine Improved_Euler

输出结果：


Improved Euler 方法

图例：


不同h每步迭代结果

Four Order Runge-Kutta Method

流程图：


原理：
![][5]
[5]: http://latex.codecogs.com/gif.latex?\begin{cases}%20y_{n+1}=y_n+\dfrac{h}{6}(K_1+2K_2+2K_3+K_4)%20\%20K_1=f(x_n,y_n)%20\%20K_2=f(x_n+\dfrac{h}{2},y_n+\dfrac{h}{2}K_1)%20\%20K_3=f(x_n+\dfrac{h}{2},y_n+\dfrac{h}{2}K_2)%20\%20K_4=f(x_n+h,y_n+hK_3)%20\end{cases}

代码：

    subroutine Four_Order_Runge_Kutta()
        implicit none
        real :: k1,k2,k3,k4
        do i=1,n
            k1=f(x(i-1),y(i-1))
            k2=f(x(i-1)+h/2,y(i-1)+h/2*k1)
            k3=f(x(i-1)+h/2,y(i-1)+h/2*k2)
            k4=f(x(i-1)+h,y(i-1)+h*k3)
            y(i)=y(i-1)+h/6*(k1+2*k2+2*k3+k4)
        end do
    end subroutine Four_Order_Runge_Kutta

输出结果：


4 Order Runge-Kutta 方法

图例：


不同h每步迭代结果

三种方法结果对比

h=0.1时
h=0.1/8时

可以看出，当h较大时，三种方法的差别还是很大的，当h逐渐减小时，三种方法的结果已基本相同。

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

推荐阅读更多精彩内容

  • LeetCode/LintCode ReviewPage 题解-总结
    背景 一年多以前我在知乎上答了有关LeetCode的问题, 分享了一些自己做题目的经验。 张土汪：刷leetcod...
    土汪阅读 12,767评论 0赞 33
  • 无标题文章
    转至元数据结尾创建： 董潇伟，最新修改于： 十二月 23, 2016 转至元数据起始第一章:isa和Class一....
    40c0490e5268阅读 1,757评论 0赞 9
  • Spring Cloud
    Spring Cloud为开发人员提供了快速构建分布式系统中一些常见模式的工具（例如配置管理，服务发现，断路器，智...
    卡卡罗2017阅读 134,845评论 18赞 139
  • 一切还长着呢
    午饭过后，D小姐照常跟着同事们在办公桌上聊会儿天，那是一天之中她感觉最轻松的时刻。可以听着这些长辈们怎么说世...
    美的岁月阅读 231评论 0赞 0
  • Sequelize(2) 定义表格
    定义 定义模型与表之间的映射,使用 define 方法. Sequelize 会自动增加 createdAt 和 ...
    kayorl阅读 7,980评论 2赞 4