回 帖 发 新 帖 刷新版面

主题:[讨论]哪位牛人帮忙本人菜鸟 关于常微分方程的!!!!

不知道如何处理这个常微分方程组 高手指点 十分感激

回复列表 (共1个回复)

沙发

给你个例子
!
!        四阶龙格库塔法求解常微分方程组
!        dx/dt = y
!        dy/dt = -x
!        x(0) = 0
!        y(0) = 1
!        其精确解为 x=sin(t),y=cos(t)
!        
!
!-----------------------------------------------------------------------
!      1. 主程序
!-----------------------------------------------------------------------
      program main
         implicit none
         real(4) :: h,t,x,y
         integer(4) :: n
         h = .05
         t = 0.
         x = 0.
         y = 1.
         do n=1,100
            call runge_kutta(h,t,x,y)
            t = n*h
            if( mod(n,10).eq.0 ) write(*,*) t,x,sin(t)
         enddo
         
      end program main

!-----------------------------------------------------------------------
!      2. 龙格库塔子程序
!-----------------------------------------------------------------------
      subroutine runge_kutta(h,t,x,y)
         implicit none
         real(4) :: h,t,x,y
         real,external :: f,g
         real(4) :: FF(4),GG(4)
         
         FF(1) = f(t,x,y)
         GG(1) = g(t,x,y)
         FF(2) = f(t+.5*h,x+.5*h*FF(1),y+.5*h*GG(1))
         GG(2) = g(t+.5*h,x+.5*h*FF(1),y+.5*h*GG(1))
         FF(3) = f(t+.5*h,x+.5*h*FF(2),y+.5*h*GG(2))
         GG(3) = g(t+.5*h,x+.5*h*FF(2),y+.5*h*GG(2))
         FF(4) = f(t+h,x+h*FF(3),y+h*GG(3))
         GG(4) = g(t+h,x+h*FF(3),y+h*GG(3))
         x = x+h/6.*(FF(1)+2.*FF(2)+2.*FF(3)+FF(4))
         y = y+h/6.*(GG(1)+2.*GG(2)+2.*GG(3)+GG(4))
         
         return
      end subroutine runge_kutta
!-----------------------------------------------------------------------
!      3. 微分函数dx/dt
!-----------------------------------------------------------------------
      function f(a,b,c)
         implicit none
         real(4) :: f,a,b,c
         f = c
      end function f
!-----------------------------------------------------------------------
!      4. 微分函数dy/dt
!-----------------------------------------------------------------------
      function g(a,b,c)
         implicit none
         real(4) :: g,a,b,c
         g = -b
      end function g

我来回复

您尚未登录,请登录后再回复。点此登录或注册