Fortran微分方程

Question

今天，我试图对Fortran 90中气体中内部能量的微分方程进行评估：

du / dt = dT / dt =-λ/ρ

其中u是内能，λ是冷却函数(它们都是温度T的函数)。ρ是质量密度，我们可以假定它是常数。

我使用的是Runge二阶方法(heun)，我确信我正确地编写了实际的求解算法，但我很确定我正在搞砸实现。我也不知道如何有效地选择任意的能量尺度。

我正在用这个子例程实现右侧：

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

考虑到lambda和rho都是n维数组，我将结果保存在一个名为u_tilde的数组中，在T= 1,000,000 K处选择一个起始条件

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

这给了我一个奇怪的温度随时间变化的图表。

​

​

我想要启动温度为一百万开尔文，然后让它冷却到10.000 K，看看它需要多长时间。如何落实这些边界条件？在RHS和在程序中设置计算循环时，我做错了什么？

Answer 1

dydx的实现只分配第一个元素。

此外，不需要为每个步骤定义循环，因为Fortran90可以执行向量操作。

对于模块化设计，我建议实现一个保存模型数据的自定义类型，比如质量密度和冷却系数。

下面是一个简单的实现示例，它只保存一个标量值，例如y‘= -c y

module mod_diffeq
use, intrinsic :: iso_fortran_env, wp => real64
implicit none

type :: model
    real(wp) :: coefficient
end type

contains

pure function dxdy(arg, x, y) result(yp)
type(model), intent(in) :: arg
real(wp), intent(in) :: x, y(:)
real(wp) :: yp(size(y))    
    yp = -arg%coefficient*y    
end function

pure function heun(arg, x0, y0, h) result(y)
type(model), intent(in) :: arg
real(wp), intent(in) :: x0, y0(:), h
real(wp) :: y(size(y0)), k0(size(y0)), k1(size(y0))
    k0 = dxdy(arg, x0, y0)
    k1 = dxdy(arg, x0+h, y0 + h*k0)        
    y = y0 + h*(k0+k1)/2
end function
end module

并利用上述模块进行了一些冷却模拟。

program FortranCoolingConsole1
use mod_diffeq
implicit none

integer, parameter :: neq = 100
integer, parameter :: nsteps = 256
! Variables
type(model):: gas
real(wp) :: x, y(neq), x_end, h
integer :: i

! Body of Console1
gas%coefficient = 1.0_wp
x = 0.0_wp
x_end = 10.0_wp
do i=1, neq
    if(i==1) then
        y(i)  = 1000.0_wp
    else
        y(i) = 0.0_wp
    end if
end do
print '(1x," ",a22," ",a22)', 'x', 'y(1)'
print '(1x," ",g22.15," ",g22.15)', x, y(1)
! Initial Conditions
h = (x_end  - x)/nsteps
! Simulation
do while(x<x_end)
    x = x + h
    y = heun(gas, x, y, h)
    print '(1x," ",g22.15," ",g22.15)', x, y(1)
end do

end program 

注意，我只跟踪neq组件的y的第一个元素。

样本输出从1000开始呈指数衰减。

                       x                   y(1)
    0.00000000000000       1000.00000000000
   0.390625000000000E-01   961.700439453125
   0.781250000000000E-01   924.867735244334
   0.117187500000000       889.445707420492
   0.156250000000000       855.380327695983
   0.195312500000000       822.619637044785
   0.234375000000000       791.113666448740
   0.273437500000000       760.814360681126
   0.312500000000000       731.675505009287
   0.351562500000000       703.652654704519
   0.390625000000000       676.703067251694
   0.429687500000000       650.785637155231
   0.468750000000000       625.860833241968
   0.507812500000000       601.890638365300
   0.546875000000000       578.838491418631
   0.585937500000000       556.669231569681
   ...

此外，如果您希望上面的代码实现runge-kutta第4顺序，您可以在mod_diffeq模块中包含以下内容

pure function rk4(arg, x0, y0, h) result(y)
type(model), intent(in) :: arg
real(wp), intent(in) :: x0, y0(:), h
real(wp) :: y(size(y0)), k0(size(y0)), k1(size(y0)), k2(size(y0)), k3(size(y0))
    k0 = dxdy(arg, x0, y0)
    k1 = dxdy(arg, x0+h/2, y0 + (h/2)*k0)
    k2 = dxdy(arg, x0+h/2, y0 + (h/2)*k1)
    k3 = dxdy(arg, x0+h, y0 + h*k2)
    y = y0 + (h/6)*(k0+2*k1+2*k2+k3)
end function