主题:[讨论]哪位牛人帮忙本人菜鸟 关于常微分方程的!!!!
伊风-私
[专家分:0] 发布于 2010-06-28 10:15:00
不知道如何处理这个常微分方程组 高手指点 十分感激
沙发
vincixy [专家分:0] 发布于 2010-08-08 22:59:00
给你个例子
!
! 四阶龙格库塔法求解常微分方程组
! 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