一个算微分方程程序求助

yifanxing1992 · 发表于 2016-7-8 13:10:13

module m_rk4
!四阶龙格库塔方法计算微分方程
contains
subroutine rk4(t0,tt,y0,N)
implicit real*8(a-z)
integer::N,i
!步长
H=(tt-t0)/N
t=t0
y=y0
do i=1,N
      call fun1(k1,t,y)
      call fun1(k2,t+H/2,y+H/2*k1)
      call fun1(k3,t+H/2,y+H/2*k2)
      call fun1(k4,t+H,y+H*k3)
      y=y+(k1+2*k2+2*k3+k4)*H/6
      t=t0+i*H
      write(11,21)t,y
end do
21  format (f15.5,2f20.8)
end subroutine rk4

subroutine fun1(f,t,y)
!待计算的函数模块，其中f(t,y)
implicit real*8(a-z)
real*8::pi
real*8::S1,S2,p1,V1,l,h,Tk1,Tk2,T1,T2,rou,g
pi=3.141592653

S1=0.36*pi*10**(-6)
S2=2.31*pi*10**(-5)
p1=1.01*10**5
l=0.345
h=0.15
V1=pi*(4.81*10**(-3))**2*(l-h)
T1=20
T2=80
Tk1=T1+273.15
Tk2=T2+273.15
rou=10**3
g=9.8

f=rou**-1*((S1/S2)-1)**-1*((p1*V1*Tk2*S1)/(Tk1*(V1+S1*y))+rou*g*(1-S1/S2))
! f=sqrt(rou**-1*((S1/S2)**2-1)**-1*(2*(p1-(p1*V1*Tk2)/(Tk1*(V1+S1*y)))+2*rou*g*(h1+(1-S1/S2)*y)))
end subroutine fun1
end module m_rk4

program rspq
!龙格库塔方法主程序

use m_rk4
implicit real*8(a-z)
integer::N=1000
open(unit=11,file='result.txt')
t0=0
t1=1
y0=0
write(11,66)
66 format(/,T5,'计算结果为:',/)
call rk4(t0,t1,y0,N)

end program rspq

以上是原程序，初始条件和方程验证无书写错误，请问为何最后的数值结果一直显示y=NAN？？？？

li913 · 发表于 2016-7-8 17:01:03

NAN = not a number.
计算过程中出现了数学错误，比如，我是说比如，对负数开方。
修改设置，使之报错，这样就知道是哪里出问题了。

[微分方程] 一个算微分方程程序求助