Fortran 90求解气体内能微分方程的实现问题咨询

我正在尝试用Fortran 90求解气体内能相关的微分方程：
$$\frac{du}{dt} = \frac{dT}{dt} = -\frac{\lambda}{\rho}$$
其中$u$为内能，$\lambda$为冷却函数（二者均仅为温度$T$的函数），$\rho$为质量密度且可视为恒定值。我采用二阶Runge-Kutta法（Heun法）求解，确认算法逻辑正确，但怀疑实现存在问题，同时不清楚如何高效选择任意能量尺度。

已实现代码 #

右端项（RHS）子例程

MODULE RHS
! right hand side
  IMPLICIT NONE
  CONTAINS
  
  SUBROUTINE dydx(neq, y, f)
    INTEGER, INTENT(IN) :: neq
    REAL*8, DIMENSION(neq), INTENT(IN) :: y
    REAL*8, DIMENSION(neq), INTENT(OUT) :: f
    
    f(1) = -y(1)

  END SUBROUTINE dydx
  
END MODULE RHS

Heun算法实现

SUBROUTINE heun(neq, h, yold, ynew)
    INTEGER, INTENT(IN) :: neq
    REAL*8, INTENT(IN) :: h 
    REAL*8, DIMENSION(neq), INTENT(IN) ::yold
    REAL*8, DIMENSION(neq), INTENT(OUT) :: ynew

    REAL*8, DIMENSION(neq) :: f, ftilde
    INTEGER :: i

    CALL dydx(neq, yold, f)

    DO i=1, neq
      ynew(i) = yold(i) + h*f(i)
    END DO

    CALL dydx(neq, ynew, ftilde)

    DO i=1, neq
      ynew(i) = yold(i) + 0.5d0*h*(f(i) + ftilde(i))
    END DO
    
  END SUBROUTINE heun

主程序部分

h = 1.d0/n
u_tilde(1) = lambda(n)/density(n)   ! lambda(3) is at about T=one million

DO i = 2, n
CALL heun(1, h*i, u_tilde(i-1), u_tilde(i))
ENDDO

运行后得到异常的温度-时间曲线。我需要设置初始温度为1,000,000K、冷却至10,000K的边界条件，想咨询当前实现中的错误点，以及如何正确设置这些边界条件。

1. 右端项子例程完全偏离物理方程 #

你的微分方程是$\frac{dT}{dt} = -\frac{\lambda(T)}{\rho}$，但当前dydx子例程返回f(1) = -y(1)，对应$\frac{dT}{dt} = -T$，和实际物理模型完全不符。

修正方法 #

需要根据当前温度值（y(1)）匹配对应的冷却函数$\lambda$，再计算正确的右端项。如果$\lambda$是离散数组（对应不同温度点），需补充插值逻辑：

MODULE RHS
  IMPLICIT NONE
  CONTAINS
  
  SUBROUTINE dydx(neq, y, f, lambda_arr, temp_arr, rho, n_lambda)
    INTEGER, INTENT(IN) :: neq, n_lambda
    REAL*8, DIMENSION(neq), INTENT(IN) :: y
    REAL*8, DIMENSION(n_lambda), INTENT(IN) :: lambda_arr, temp_arr
    REAL*8, INTENT(IN) :: rho
    REAL*8, DIMENSION(neq), INTENT(OUT) :: f
    
    REAL*8 :: lambda_current
    ! 调用插值子例程获取当前温度对应的lambda值
    CALL linear_interp(y(1), temp_arr, lambda_arr, n_lambda, lambda_current)
    
    f(1) = -lambda_current / rho
  END SUBROUTINE dydx

  ! 线性插值子例程
  SUBROUTINE linear_interp(x, x_arr, y_arr, n, y_out)
    INTEGER, INTENT(IN) :: n
    REAL*8, INTENT(IN) :: x
    REAL*8, DIMENSION(n), INTENT(IN) :: x_arr, y_arr
    REAL*8, INTENT(OUT) :: y_out
    
    INTEGER :: i
    IF(x <= x_arr(1)) THEN
      y_out = y_arr(1)
    ELSE IF(x >= x_arr(n)) THEN
      y_out = y_arr(n)
    ELSE
      DO i=1, n-1
        IF(x >= x_arr(i) .AND. x <= x_arr(i+1)) THEN
          y_out = y_arr(i) + (x - x_arr(i))*(y_arr(i+1)-y_arr(i))/(x_arr(i+1)-x_arr(i))
          EXIT
        END IF
      END DO
    END IF
  END SUBROUTINE linear_interp
  
END MODULE RHS

2. 时间步长调用逻辑错误 #

主程序中调用heun时传入h*i，导致步长随循环次数线性递增（第2步用2h，第3步用3h...），完全违背数值积分的步长规则——Heun法每一步应使用固定小步长（或自适应步长）。

修正方法 #

定义总模拟时间，计算固定步长，循环中每次传入相同步长：

REAL*8, DIMENSION(n) :: t_arr, T_arr
REAL*8 :: h, t_total
INTEGER :: i

t_total = 1.0d6 ! 根据物理需求设置总时间
h = t_total / n
T_arr(1) = 1.0d6 ! 初始温度1e6 K
t_arr(1) = 0.0d0

DO i = 2, n
  CALL heun(1, h, T_arr(i-1), T_arr(i))
  t_arr(i) = t_arr(i-1) + h
ENDDO

3. 初始条件设置错误 #

当前代码中u_tilde(1) = lambda(n)/density(n)完全不符合初始温度要求，初始条件应直接赋值温度（或对应内能）。

修正方法 #

直接将初始温度赋值给数组首元素：

T_arr(1) = 1.0d6 ! 初始温度1,000,000 K

若模拟内能，可根据理想气体关系$u = C_v T$（$C_v$为定容比热）计算初始内能。

4. 冷却终止条件设置 #

不需要固定循环次数n，应设置温度低于10,000 K时终止积分：

REAL*8 :: T_current, T_next, t_current, h
INTEGER :: step

T_current = 1.0d6
t_current = 0.0d0
h = 100.0d0 ! 初始步长，可根据精度调整

OPEN(UNIT=10, FILE="cooling_curve.dat", STATUS="REPLACE")
WRITE(10,*) t_current, T_current

DO WHILE(T_current > 1.0d4)
  CALL heun(1, h, T_current, T_next)
  ! 可选：自适应调整步长，保证计算精度
  IF(ABS(T_next - T_current)/T_current > 1.0d-3) THEN
    h = h * 0.5d0
    CYCLE
  END IF
  t_current = t_current + h
  T_current = T_next
  WRITE(10,*) t_current, T_current
END DO
CLOSE(10)

5. 能量尺度选择方案 #

通过无量纲化处理实现任意能量尺度适配：

• 选择特征温度$T_0$（比如初始温度1e6 K），特征时间$\tau = \rho T_0 / \lambda_0$（$\lambda_0$为$T_0$对应的冷却函数值）。
• 定义无量纲温度$\theta = T/T_0$，无量纲时间$\tau' = t/\tau$，方程转化为：
  $$\frac{d\theta}{d\tau'} = -\frac{\lambda(\theta T_0)}{\lambda_0}$$
  该处理可将变量缩放到0~1范围，避免数值溢出，同时方便调整尺度。

内容的提问来源于stack exchange，提问作者Marco Leonardi