使用Runge-Kutta求解Bernoulli方程的Fortran代码报错排查

原代码与编译报错 #

原Fortran代码 #

program Bernoulli_Equation
    integer N ,i,u(i+1),p(i+1),rho(i+1),u(i),p(i),rho(i),u(N), p(N), rho(N) ! Nombre d'itérations
    real t , k1,l1,m1,k2,l2,m2,k3,l3,m3,k4,l4,m4 ,dt! Variables dépendantes
    
    ! Conditions initiales
    u(1) = 0.0
    p(1) = 1.0
    rho(1) = 1.0
    N=100
    dt=0.01
    open(unit=5,file='OHRK4.txt')
    ! Boucle de calcul
    do i = 1, N-1
        t = i * dt
        
        ! Calcul des pentes
        k1 = dt * du_dt(u(i), p(i), rho(i))
        l1 = dt * dp_dt(u(i), p(i), rho(i))
        m1 = dt * drho_dt(u(i), p(i), rho(i))
        
        k2 = dt * du_dt(u(i) + 0.5*k1, p(i) + 0.5*l1, rho(i) + 0.5*m1)
        l2 = dt * dp_dt(u(i) + 0.5*k1, p(i) + 0.5*l1, rho(i) + 0.5*m1)
        m2 = dt * drho_dt(u(i) + 0.5*k1, p(i) + 0.5*l1, rho(i) + 0.5*m1)
        
        k3 = dt * du_dt(u(i) + 0.5*k2, p(i) + 0.5*l2, rho(i) + 0.5*m2)
        l3 = dt * dp_dt(u(i) + 0.5*k2, p(i) + 0.5*l2, rho(i) + 0.5*m2)
        m3 = dt * drho_dt(u(i) + 0.5*k2, p(i) + 0.5*l2, rho(i) + 0.5*m2)
        
        k4 = dt * du_dt(u(i) + k3, p(i) + l3, rho(i) + m3)
        l4 = dt * dp_dt(u(i) + k3, p(i) + l3, rho(i) + m3)
        m4 = dt * drho_dt(u(i) + k3, p(i) + l3, rho(i) + m3)
        
        ! Mise à jour des variables dépendantes
        u(i+1) = u(i) + (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0
        p(i+1) = p(i) + (l1 + 2.0*l2 + 2.0*l3 + l4) / 6.0
        rho(i+1) = rho(i) + (m1 + 2.0*m2 + 2.0*m3 + m4) / 6.0
    end do
    
    ! Affichage des résultats
    do i = 1, N
        t = i * dt
        write(*,*) t, u(i), p(i), rho(i)
    end do
    
contains
    
    ! Fonctions pour calculer les

编译报错信息 #

2 |     integer N ,i,u(i+1),p(i+1),rho(i+1),u(i),p(i),rho(i),u(N), p(N), rho(N) ! Nombre d'itérations
      |                        1
Error: Explicit shaped array with nonconstant bounds at (1)
bernouli.f90:6:4:

    6 |     u(1) = 0.0
      |    1
Error: The function result on the lhs of the assignment at (1) must have the pointer attribute.
bernouli.f90:7:4:

    7 |     p(1) = 1.0
      |    1
Error: The function result on the lhs of the assignment at (1) must have the pointer attribute.
bernouli.f90:8:4:

    8 |     rho(1) = 1.0
      |    1
Error: The function result on the lhs of the assignment at (1) must have the pointer attribute.
bernouli.f90:34:8:

   34 |         u(i+1) = u(i) + (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0
      |        1
Error: The function result on the lhs of the assignment at (1) must have the pointer attribute.
bernouli.f90:35:8:

   35 |         p(i+1) = p(i) + (l1 + 2.0*l2 + 2.0*l3 + l4) / 6.0
      |        1
Error: The function result on the lhs of the assignment at (1) must have the pointer attribute.
bernouli.f90:36:8:

   36 |         rho(i+1) = rho(i) + (m1 + 2.0*m2 + 2.0*m3 + m4) / 6.0
      |        1
Error: The function result on the lhs of the assignment at (1) must have the pointer attribute.

错误原因分析 #

1. 数组声明严重错误：
   • 同一行反复定义u、p、rho为不同大小的数组，Fortran不允许重复定义变量。
   • 使用未初始化的i、N作为数组边界，显式形状数组要求边界必须是编译期常量，不能用运行时才赋值的变量。
   • 错误地将实数数组声明为integer类型，导致后续赋值实数时类型不匹配。
2. 数组初始化顺序错误：先给数组元素赋值，再给N赋值，此时数组还未正确分配空间。
3. 导数函数未实现：contains块内的du_dt、dp_dt、drho_dt函数只有注释，没有具体实现，编译时会找不到函数定义。

修正后的代码 #

program Bernoulli_Equation
    implicit none  ! 强制显式声明所有变量，避免隐式类型错误
    integer :: N, i
    real :: t, dt, k1, l1, m1, k2, l2, m2, k3, l3, m3, k4, l4, m4
    real, allocatable :: u(:), p(:), rho(:)  ! 声明可分配数组
    
    ! 先设置迭代次数与时间步长
    N = 100
    dt = 0.01
    
    ! 分配数组空间
    allocate(u(N), p(N), rho(N))
    
    ! 初始条件
    u(1) = 0.0
    p(1) = 1.0
    rho(1) = 1.0
    
    open(unit=5, file='OHRK4.txt')
    
    ! Runge-Kutta迭代计算
    do i = 1, N-1
        t = i * dt
        
        ! 计算各阶斜率
        k1 = dt * du_dt(u(i), p(i), rho(i))
        l1 = dt * dp_dt(u(i), p(i), rho(i))
        m1 = dt * drho_dt(u(i), p(i), rho(i))
        
        k2 = dt * du_dt(u(i) + 0.5*k1, p(i) + 0.5*l1, rho(i) + 0.5*m1)
        l2 = dt * dp_dt(u(i) + 0.5*k1, p(i) + 0.5*l1, rho(i) + 0.5*m1)
        m2 = dt * drho_dt(u(i) + 0.5*k1, p(i) + 0.5*l1, rho(i) + 0.5*m1)
        
        k3 = dt * du_dt(u(i) + 0.5*k2, p(i) + 0.5*l2, rho(i) + 0.5*m2)
        l3 = dt * dp_dt(u(i) + 0.5*k2, p(i) + 0.5*l2, rho(i) + 0.5*m2)
        m3 = dt * drho_dt(u(i) + 0.5*k2, p(i) + 0.5*l2, rho(i) + 0.5*m2)
        
        k4 = dt * du_dt(u(i) + k3, p(i) + l3, rho(i) + m3)
        l4 = dt * dp_dt(u(i) + k3, p(i) + l3, rho(i) + m3)
        m4 = dt * drho_dt(u(i) + k3, p(i) + l3, rho(i) + m3)
        
        ! 更新变量
        u(i+1) = u(i) + (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0
        p(i+1) = p(i) + (l1 + 2.0*l2 + 2.0*l3 + l4) / 6.0
        rho(i+1) = rho(i) + (m1 + 2.0*m2 + 2.0*m3 + m4) / 6.0
        
        ! 将结果写入文件
        write(5,*) t, u(i), p(i), rho(i)
    end do
    
    ! 输出最后一个时间步的结果
    write(5,*) N*dt, u(N), p(N), rho(N)
    
    ! 屏幕输出结果
    do i = 1, N
        t = i * dt
        write(*,*) t, u(i), p(i), rho(i)
    end do
    
    ! 释放数组空间
    deallocate(u, p, rho)
    close(5)
    
contains
    ! 请根据你的Bernoulli方程实际形式修改以下导数函数
    real function du_dt(u_val, p_val, rho_val)
        real, intent(in) :: u_val, p_val, rho_val
        ! 示例：假设du/dt = -p_val/rho_val （替换为你的实际方程）
        du_dt = -p_val / rho_val
    end function du_dt
    
    real function dp_dt(u_val, p_val, rho_val)
        real, intent(in) :: u_val, p_val, rho_val
        ! 示例：假设dp/dt = -rho_val * u_val （替换为你的实际方程）
        dp_dt = -rho_val * u_val
    end function dp_dt
    
    real function drho_dt(u_val, p_val, rho_val)
        real, intent(in) :: u_val, p_val, rho_val
        ! 示例：假设drho/dt = 0 （替换为你的实际方程）
        drho_dt = 0.0
    end function drho_dt

end program Bernoulli_Equation

关键修正点说明 #

• 添加implicit none：禁止隐式变量声明，避免类型错误。
• 修正数组声明：使用可分配数组real, allocatable :: u(:), p(:), rho(:)，先给N赋值后再用allocate分配空间，符合Fortran数组声明规则。
• 修正变量类型：将u、p、rho改为real类型，与后续赋值的实数匹配。
• 补全导数函数：在contains块内实现了三个导数函数的示例，需根据你实际的Bernoulli方程形式修改函数体。
• 调整执行顺序：先设置N和dt，再分配数组，最后赋值初始条件，避免数组未分配就使用的错误。
• 添加数组释放与文件关闭：养成良好的内存管理和文件操作习惯。

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