你的提示是c++出问题，不合常理，可能是软件问题。我这里是数组越界，

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

den(-3:103,0:W-1,3)
。。。
do while(n<=W-1)
。。。
 den(j,n+1,3)=1d0/3d0。。。

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

program main
    implicit none
    integer i,j,W,n
    real(8) gamma,ux,delt_t,delt_x,total_t,uxp,uxm
    real(8),allocatable::den(:,:,:),u(:,:,:),p(:,:,:),den_u(:,:,:),e(:,:,:)
    real(8) f1_p(-3:103),f2_p(-3:103),f3_p(-3:103),f1_m(-3:103),f2_m(-3:103),f3_m(-3:103)
    gamma=1.4d0
    total_t=0.14d0
    delt_x=0.01d0
    delt_t=0.001d0
    W=total_t/delt_t+1
    allocate(den(-3:103,0:W-1,3),u(-3:103,0:W-1,3),p(-3:103,0:W-1,3),den_u(-3:103,0:W-1,3),e(-3:103,0:W-1,3))

    n=0
    !----initial value
    do j=-3,103
        u(j,n,3)=0d0
        if(j*delt_x<0.5d0)then
            den(j,n,3)=1d0
            p(j,n,3)=1d0
            else
                den(j,n,3)=0.125d0
                p(j,n,3)=0.1d0
        end if
        den_u(j,n,3)=den(j,n,3)*u(j,n,3)
        e(j,n,3)=p(j,n,3)/(gamma-1d0)+0.5d0*den(j,n,3)*u(j,n,3)**2
    end do

    !----runge-kutta
    do while(n<=W-1)

        !----runge-kutta step1
        i=3
        call FVS(i,u,p,den,n,f1_p,f1_m,f2_p,f2_m,f3_p,f3_m,W)

        do j=0,100
            ux=uxp(f1_p,j)+uxm(f1_m,j)
            den(j,n,1)=den(j,n,3)-delt_t*ux

            ux=uxp(f2_p,j)+uxm(f2_m,j)
            den_u(j,n,1)=den_u(j,n,3)-delt_t*ux
            u(j,n,1)=den_u(j,n,1)/den(j,n,1)

            ux=uxp(f3_p,j)+uxm(f3_m,j)
            e(j,n,1)=e(j,n,3)-delt_t*ux
            p(j,n,1)=(gamma-1d0)*(e(j,n,1)-0.5d0*den(j,n,1)*u(j,n,1)**2)
        end do

        do j=-3,-1
            u(j,n,1)=u(0,n,1)
            den(j,n,1)=den(0,n,1)
            p(j,n,1)=p(0,n,1)
        end do
        do j=101,103
            u(j,n,1)=u(100,n,1)
            den(j,n,1)=den(100,n,1)
            p(j,n,1)=p(100,n,1)
        end do

        !----runge-kutta step2
        i=1
        call FVS(i,u,p,den,n,f1_p,f1_m,f2_p,f2_m,f3_p,f3_m,W)

        do j=0,100
            ux=uxp(f1_p,j)+uxm(f1_m,j)
            den(j,n,2)=3d0/4d0*den(j,n,3)+1d0/4d0*(den(j,n,1)-delt_t*ux)

            ux=uxp(f2_p,j)+uxm(f2_m,j)
            den_u(j,n,2)=3d0/4d0*den_u(j,n,3)+1d0/4d0*(den_u(j,n,1)-delt_t*ux)
            u(j,n,2)=den_u(j,n,2)/den(j,n,2)

            ux=uxp(f3_p,j)+uxm(f3_m,j)
            e(j,n,2)=3d0/4d0*e(j,n,3)+1d0/4d0*(e(j,n,1)-delt_t*ux)
            p(j,n,2)=(gamma-1d0)*(e(j,n,2)-0.5d0*den(j,n,2)*u(j,n,2)**2)
        end do

        do j=-3,-1
            u(j,n,2)=u(0,n,2)
            den(j,n,2)=den(0,n,2)
            p(j,n,2)=p(0,n,2)
        end do
        do j=101,103
            u(j,n,2)=u(100,n,2)
            den(j,n,2)=den(100,n,2)
            p(j,n,2)=p(100,n,2)
        end do

        !----runge-kutta step3
        i=2
        call FVS(i,u,p,den,n,f1_p,f1_m,f2_p,f2_m,f3_p,f3_m,W)

        do j=0,100
            ux=uxp(f1_p,j)+uxm(f1_m,j)
            den(j,n+1,3)=1d0/3d0*den(j,n,3)+2d0/3d0*(den(j,n,2)-delt_t*ux)

            ux=uxp(f2_p,j)+uxm(f2_m,j)
            den_u(j,n+1,3)=1d0/3d0*den_u(j,n,3)+2d0/3d0*(den_u(j,n,2)-delt_t*ux)
            u(j,n+1,3)=den_u(j,n+1,3)/den(j,n+1,3)

            ux=uxp(f3_p,j)+uxm(f3_m,j)
            e(j,n+1,3)=1d0/3d0*e(j,n,3)+2d0/3d0*(e(j,n,2)-delt_t*ux)
            p(j,n+1,3)=(gamma-1d0)*(e(j,n+1,3)-0.5d0*den(j,n+1,3)*u(j,n+1,3)**2)
        end do

        do j=-3,-1
            u(j,n+1,3)=u(0,n+1,3)
            den(j,n+1,3)=den(0,n+1,3)
            p(j,n+1,3)=p(0,n+1,3)
        end do
        do j=101,103
            u(j,n+1,3)=u(100,n+1,3)
            den(j,n+1,3)=den(100,n+1,3)
            p(j,n+1,3)=p(100,n+1,3)
        end do

    n=n+1
    end do

    open(unit=10,file="FVS.txt")
    do j=0,100
        write(10,"(4E15.7)") j*delt_x,den(j,W-1,3),u(j,W-1,3),p(j,W-1,3)
    end do
    close(unit=10)

end

subroutine FVS(i,u,p,den,n,f1_p,f1_m,f2_p,f2_m,f3_p,f3_m,W)
    implicit none
    integer i,j,W,n
    real(8) gamma,epsilon,wp,wm,coe
    real(8) c,lamda1,lamda1_p,lamda1_m,lamda2,lamda2_p,lamda2_m,lamda3,lamda3_p,lamda3_m
    real(8) f1_p(-3:103),f2_p(-3:103),f3_p(-3:103),f1_m(-3:103),f2_m(-3:103),f3_m(-3:103)
    real(8) u(-3:103,0:W-1,3),den(-3:103,0:W-1,3),p(-3:103,0:W-1,3)
    gamma=1.4d0
    epsilon=0.1d0
    do j=-3,103
        c=sqrt(gamma*p(j,n,i)/den(j,n,i))
        lamda1=u(j,n,i)
        lamda2=u(j,n,i)-c
        lamda3=u(j,n,i)+c
        lamda1_p=(lamda1+sqrt(lamda1**2+epsilon**2))/2d0
        lamda1_m=(lamda1-sqrt(lamda1**2+epsilon**2))/2d0
        lamda2_p=(lamda2+sqrt(lamda2**2+epsilon**2))/2d0
        lamda2_m=(lamda2-sqrt(lamda2**2+epsilon**2))/2d0
        lamda3_p=(lamda3+sqrt(lamda3**2+epsilon**2))/2d0
        lamda3_m=(lamda3-sqrt(lamda3**2+epsilon**2))/2d0
        wp=(3d0-gamma)*(lamda2_p+lamda3_p)*c**2/(2d0*(gamma-1d0))
        wm=(3d0-gamma)*(lamda2_m+lamda3_m)*c**2/(2d0*(gamma-1d0))
        coe=den(j,n,i)/(2d0*gamma)
        f1_p(j)=coe*(2d0*(gamma-1d0)*lamda1_p+lamda2_p+lamda3_p)
        f1_m(j)=coe*(2d0*(gamma-1d0)*lamda1_m+lamda2_m+lamda3_m)
        f2_p(j)=coe*(2d0*(gamma-1d0)*lamda1_p*lamda1+lamda2_p*lamda2+lamda3_p*lamda3)
        f2_m(j)=coe*(2d0*(gamma-1d0)*lamda1_m*lamda1+lamda2_m*lamda2+lamda3_m*lamda3)
        f3_p(j)=coe*((gamma-1d0)*lamda1_p*lamda1**2+0.5d0*lamda2_p*lamda2**2+0.5d0*lamda3_p*lamda3**2+wp)
        f3_m(j)=coe*((gamma-1d0)*lamda1_m*lamda1**2+0.5d0*lamda2_m*lamda2**2+0.5d0*lamda3_m*lamda3**2+wm)
    end do
end subroutine

function uxp(up,j)
    implicit none
    integer j
    real(8) up(-3:103)
    real(8) uxp,delt_x
    delt_x=0.01d0
    uxp=(-2d0*up(j-3))+15d0*up(j-2)-60d0*up(j-1)+20d0*up(j)+30d0*up(j+1)-3d0*up(j+2)/(60d0*delt_x)
end function

function uxm(um,j)
    implicit none
    integer j
    real(8) um(-3:103)
    real(8) uxm,delt_x
    delt_x=0.01d0
    uxm=(3d0*um(j-2)-30d0*um(j-1)-20d0*um(j)+60d0*um(j+1)-15d0*um(j+2)+2d0*um(j+3))/(60d0*delt_x)
end function

是个求解sod激波管的问题，用了流通矢量分裂，时间推进是3阶runge-kutta，两个附件是一模一样的代码。