程序员数学之美

5. Ordinary Differential Equatio

2015-05-14  本文已影响1016人  Uri

求解问题

Write a program to solve the following ordinary differential equation by

  • basic Euler method
  • improved Euler method
  • four-order Runge-Kutta 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逐渐减小时,三种方法的结果已基本相同。

上一篇下一篇

猜你喜欢

热点阅读