Fortran Coder

查看: 4190|回复: 1

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

[复制链接]

26

帖子

8

主题

0

精华

熟手

F 币
123 元
贡献
78 点
发表于 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????

790

帖子

2

主题

0

精华

大宗师

F 币
3765 元
贡献
2255 点
发表于 2016-7-8 17:01:03 | 显示全部楼层
NAN = not a number.
计算过程中出现了数学错误,比如,我是说比如,对负数开方。
修改设置,使之报错,这样就知道是哪里出问题了。
QQ截图20160708170019.jpg
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-3-29 12:52

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表