求助四阶龙格库塔

观星游月

用的是vs，代码如下，用四阶龙格库塔和差值计算TOV方程的几组数值解，运行窗口是空白的，不知道怎么改

[Fortran] 纯文本查看 复制代码

program gongcheng

    implicit none

    real(kind=8),parameter::G=1.3218278962236387e-42

    real(kind=8),parameter::m0=1.1158032783756145e60

    real(kind=8),parameter::R=20.0e18

    real(kind=8),parameter::pi=3.1415926535

    real(kind=8),parameter::h0=1.0e15

    real(kind=8)::p0(800) 

    real,external::f

    real,external::q

    integer::i

    real(kind=8)::rho(990),ene(990),pre(990) 

    

        real(kind=8)::x(20000),z(20000),y(20000),p(20000),m(20000),e(20000)

    double precision h,n,k1,k2,k3,k4,M1,M2,M3,M4

    

    do i=1,800

        p0(i)=i*1.d0  

        call rungekutta(x,y,z,f,q)    

    write(*,*)x(i+1),z(i+1)/m0

    end do

    end program

    

    function f(r,p,m)

    implicit none

    real(kind=8)::r,p,m,e

    real::f

    real(kind=8)::rho(990),ene(990),pre(990)

    integer::i

    real(kind=8),parameter::G=1.3218278962236387e-42

    real(kind=8),parameter::pi=3.1415926535

    open(unit=26,file='EOS.txt')

    do i=1,990

       read(26,*)rho(i),ene(i),pre(i)

    end do

    call spline(pre,ene,990,p,1,e)

    f=-(e+p)*(G*m+4*pi*G*(r**3)*p)/(r*(r-2*G*m))

    return

    end

    

    function q(r,p,m)

    implicit none

    real(kind=8)::e,p,r,m

    integer::i

    real(kind=8)::rho(990),ene(990),pre(990)

    real::q

    real(kind=8),parameter::pi=3.1415926535 

    open(unit=26,file='EOS.txt')

    do i=1,990

       read(26,*)rho(i),ene(i),pre(i)

    end do

    call spline(pre,ene,990,p,1,e)

    q=4*pi*e*r**2 

    return

    end

    

    

    subroutine rungekutta(x,y,z,f,q)

    real(kind=8)::x(20000),z(20000),y(20000),p(20000),m(20000)

    integer::i

    real,external::f

    real,external::q

    double precision h,n,k1,k2,k3,k4,M1,M2,M3,M4

    n=20.0e18/1.0e15

    h=1.0e15

    x(1)=0.0

    y(1)=p0

    z(1)=0.0

    do i=1,n-1

        y(i)=p(i)

        z(i)=m(i)

        x(i+1)=(i+1)*h

        k1=f(x(i),y(i))

        M1=q(x(i),y(i))

        k2=f(x(i)+0.5*h,y(i)+0.5*h*k1)

        M2=q(x(i)+0.5*h,y(i)+0.5*h*M1)

        k3=f(x(i)+0.5*h,y(i)+0.5*h*k2)

        M3=q(x(i)+0.5*h,y(i)+0.5*h*M2)

        k4=f(x(i)+h,y(i)+h*k3)

        M4=q(x(i)+h,y(i)+h*M3)

        y(i+1)=y(i)+h*(k1+k2+k3+k4)/6

        z(i+1)=z(i)+h*(M1+M2+M3+M4)/6 

        

        if (y(i+1)<=1.0*10**(-5))then

        exit

        end if

    end do

    end subroutine rungekutta



    

    

    subroutine spline(x,y,n,sx,m,f) !x,y:差值基点坐标，n：节点个数，sx：需要返回差值数值的x坐标，m：sx的容量，f：差值结果返回值                   

    implicit none

    integer:: i, j, k, m, n, n1, n2

    Real(kind=8)::x(n),y(n),sx,f

    Real(kind=8):: s2(n), h(n), dy(n), s(n), e(n)     !直接都给数组定义N的空间大小，方便调用 

    Real(kind=8):: z, h1, h2, h3, h4

  

     n1=n-1!间距个数

     n2=n-2!内点个数

     

     do i = 1,n1

       h(i)  = x(i+1) - x(i)!计算点间距 

     dy(i) = (y(i+1) - y(i) ) / H(i)!计算斜率sj

    end do

    

     S2(1) = 0.d0; S2(N) = 0.d0!自然边界条件，两侧二阶导数为零

     

     DO i = 2, N1!求解内结点处的二阶导数

        S2(i) = 6.d0 * ( DY(i) - DY(i-1) )

     ENDDO

    

     Z = 0.5d0 / ( H(1) + H(2) )

     S(1) = -H(2) * Z

     E(1) = S2(2) * Z

     DO i = 2, N2

        k    = i - 1

        j    = i + 1

        Z    = 1.d0 / ( 2.d0*( H(i)+H(j) ) + H(i)*S(k) )

        S(i) = -H(j) * Z

        E(i) = ( S2(j)-H(i)*E(k) ) * Z

     ENDDO

    

     S2(N1) = E(N2)

     DO i = N2, 2, -1

        k     = i - 1

        S2(i) = S(k)*S2(i+1) + E(k)

     ENDDO

    

     DO i = 1, N1

        S(i) = ( S2(i+1) - S2(i) ) / H(i)

     ENDDO

    

     !具体求结果

     i = 2

     k = 1  

     DO j = 1, M

        DO 

           IF(SX<X(1).OR.SX>X(N)) THEN

            F=0.0       

            ELSEIF ( SX > X(i) ) THEN

                k = i

                i = i + 1

            ELSE

                EXIT

            ENDIF

        ENDDO

        H1    = SX - X(k)

        H2    = SX - X(i)

        H3    = H1 * H2

        H4    = S2(k) + H1*S(k)

        Z     = ( S2(i) + S2(k) + H4 ) / 6.d0

        F  = Y(k) + H1*DY(k) + H3*Z

     ENDDO    

    end subroutine spline

观星游月

fcode 发表于 2023-2-2 08:46
你先把死循环和p0没有传入的问题修改了吧。

帮助你理解需求并完全纠正逻辑问题，需要太多的精力，我不一定 ...

感谢，我先改着

fcode

你先把死循环和p0没有传入的问题修改了吧。

帮助你理解需求并完全纠正逻辑问题，需要太多的精力，我不一定能帮助你。

你可以把上述问题修改后再看看，有问题再继续帖新的代码求助。

另外，断点单步调试会对你有很大的帮助： http://debug.w.fcode.cn  和 http://v.fcode.cn/video-debug.html

观星游月

感谢回答，过程是这样，我想求一个微分方程
f=-(e+p)*(G*m+4*pi*G*(r**3)*p)/(r*(r-2*G*m))
也就是
dp/dr=-(e+p)*(G*m+4*pi*G*(r**3)*p)/(r*(r-2*G*m))
的数值解，这里面有e,p,m,r四个变量，其他的都是常量，其中又有关于变量m的方程（叫q）
dm/dr=4*pi*e*r**2
变量e和p有一个函数关系，通过一个插值程序得出，调用txt文件中的990组pre，ene作为插值基点x，y
e和p同时又和r有一个函数关系，e(r),p（r), m也有m(r)关系
,插值得出方程之后回代到方程f中得到r,p(r)的偏微分方程，接着确定步长h，给定p的一个初值p0，把r的范围分成若干份从0到20.0e18用四阶龙格库塔循环计算p,m的值，
跳出循环的条件二选一，r的值达到20.0e18，或者p的计算值小于1.0e-5
p0取1到800重复上述过程，取循环结束的r和m/m0的800组值。

这个是论文导师给的任务，Fortran实在是刚上手，不熟练

fcode

咋说呢，我也没有头绪，发现了好多个问题。
因为不知道你的意图，也不知道怎么改。

例如：
273行开始
DO j = 1, M
        DO !//此处死循环
           IF(SX<X(1).OR.SX>X(N)) THEN
死循环了


再例如。
主程序中，p0(i)=i*1.d0
但p0并没有传入rungekutta龙格库塔子程序中。该子程序中的 p0 没有定义，也没有值。

观星游月

没啥头绪，有无大佬指点