【原】OMB | 衰减问题

1 摘要

本章使用衰减问题作为第一个示例，向读者介绍数值建模的大体轮廓。讨论了有限差分、显式和隐式方案，以及一致性、准确性、稳定性和求解效率的条件。

2 衰减问题

2.1 问题概括

衰减问题，简单来说就是物质的浓度以一定速率随时间减少的问题，整个物理过程可以用下方的方程表示：

其中 是物质的浓度，是时间，（希腊符号“kappa”）表示衰减速率，是一个常量参数。当 κ 等于 0时， 保持其初始值不变；当 不等于 0 时，因为浓度始终为正数，所以右侧为负数。因此，在任意时间点，将以与浓度本身成比例的速率随时间逐渐降低。在数学中， 该方程被称为一阶常微分方程，两边同时积分，可以求出其解析解：

                                       其中是指指数函数，表示初始浓度，当，即衰减速率确定以后，任意时间点的浓度值因初始浓度的不同而不同，因此要确定唯一解，还需要初始浓度，即，此时方程变为：

当知道时的物质浓度时，利用此方程我们就可以知道任意时刻物质的浓度，得到唯一解。相应地，叫做初始条件，像上式一样带初始条件的常微分方程问题被称为初值问题。

2.2 计算机求解问题

由于计算机不能直接求解微分方程，也不能求出其解析解，所以面对微分问题，一般是利用差商取代微商的方法，通过有限差分逼近求解，因此具有一定的误差。其中，形如的式子叫做有限差分，叫做差商，叫做微商，所以有限差分法的基本原理就是，对微分方程中的微分项进行直接差分近似，从而将微分方程转化为代数方程组求解。主要的有限差分法有：向前差分、向后差分和中心差分。

2.2.1 时间级别和时间步长

使用离散时间步长 ，我们可以将 衰减方程表述为：

其中整数 指的是某个时间级别。这个时间指数不要与幂函数的指数相混淆。通常，给出模拟开始时的浓度， 表示一个时间步长后的浓度（持续时间为 ）， 表示两个时间步长后的浓度，依此类推。

2.2.2 显式格式

将未知变量移向等式左侧，已知变量移向等式右侧， 衰减方程就变为：

其中 是指需要与 和 值一起规定的初始浓度。这种迭代方法使用在某个时间级别 已知的值来预测下一个时间级别 的 值，因此称为显式格式。

从中可以看出，对于每个时间步长，浓度都会以迭代方式减少一定比例。该部分由给出。但由于物质的浓度不能小于0，所以要求，否则预测的浓度将变为负数，这将没有意义。也就是说：

这个被称之为数值稳定条件。因此，显示格式衰减方程的预测只有在数值稳定条件满足时才是稳定的。因此，可以选择的最大时间步长取决于 的值。

2.2.3 隐式格式

若是在向后差分，则得到：

因为等式两边都有未知项，因此称为隐式格式，将已知项和未知项分离，则得到：

与显式方法相比，这种隐式方案的明显优势在于它对于任何 值都是数值稳定的。后面方程式中的分母总是大于 1，因此浓度会随时间逐渐减小（并且符号不变）。

2.2.4 混合格式

还可以混合使用显式和隐式方案，可以表述为：

其中权重因子（希腊符号“alpha”）必须从零到 1 之间的范围内选择。选择给出完全隐式方案，而导致完全显式方案。当时，我们得到一个所谓的半隐式方案。

2.2.5 其他格式

除了以上常用格式外，还有“Runge-Kutta scheme”以及“Adams-Bashforth scheme”等格式，相关内容在这没有讨论，感兴趣的话可自行探索。不同差分格式以及时间步长的选择对预测模型的准确性和效率有一定影响。

利用计算机进行差分格式求解时，除了上文显式格式特指的数值稳定条件外，还需满足一下条件：

1. 一致性条件 初始浓度的衰减问题的精确解析解由下式给出：

   如果数值模型的有限差分解在数值时间步长变得微乎其微时收敛于控制微分方程的解，则称数值模型是一致的。这意味着随着时间步长的减少，我们的模型预测的浓度应该更接近真实解。
2. 精度条件 在使用有限差分时会产生称为截断误差的特定误差。舍入误差是另一个误差来源，与计算机只能用有限位数表示数字这一事实有关。在模拟期间，这两个错误都应该保持在相当小的水平
3. 效率条件 大型程序可能需要大量计算机空间用于数据输出和存储，并且模型运行的完成可能需要很长时间。因此，必须以高效的方式编写模型代码，以便在合理的时间跨度内完成任务，并且不会用大量数据“塞满”计算机。

2.3 代码参考

前提条件：初始浓度为 100% ，，时间步长或，模拟15个小时的物质浓度变化，并探究显式格式和隐式格式的不同。语言选择：计算本部分：Fortran;  画图部分：Python
程序代码：

PROGRAM decay

!**********************************************
!* Fortran code dealing with the decay problem. 
!*
!* Author: J. Kaempf, 2008
!* Change: Daimu, 2023
!**********************************************

! Text after pronounciation marks are treated as comments.
! Every FORTRAN program starts with the "program" statement.
! The program name, "decay" here, should be different from
! the file name. Lower case or upper case doesn't matter in
! FORTRAN.
! 
! Now comes the declaration section.
! Every symbol used later on needs to be declared.

INTEGER :: n    ! time level
INTEGER :: ntot ! total number of iterations
INTEGER :: nout ! output every nout's iteration step
INTEGER :: mode ! use either explicit or implicit scheme 

REAL :: C_explicit       ! concentration at time level n
REAL :: CN_implicit      ! "N" for new, concentration at next time level n+1
REAL :: C_implicit       ! concentration at time level n
REAL :: CN_explicit      ! "N" for new, concentration at next time level n+1
REAL :: CTRUE   ! exact analytical solution for comparison
REAL :: CZERO   ! Initial value of C
REAL :: dt      ! time step
REAL :: kappa   ! decay constant
REAL :: time    ! time counter
REAL :: fac_explicit     ! factor used in the scheme
REAL :: fac_implicit     ! factor used in the scheme

! Now comes the initialization of variables and parameters.

CZERO = 100.0   ! initial concentration is 100
C_explicit = CZERO       ! initial value is used for exact solution
C_implicit = CZERO
kappa = 0.0001  ! value of decay constant; you can also write this as kappa = 1.0e-4
dt = 3600.      ! value of time step; dt = 60 seconds here

! Here is an example of an IF-statement that writes a message to the command screen,
! if the stability condition is violated.

!mode = 2        ! mode = 1 (explicit scheme); mode = 2 (implicit scheme)

fac_explicit = 1.0-dt*kappa

fac_implicit = 1.0/(1.+dt*kappa)

IF(mode == 1)THEN
   if(fac<=0.0)write(6,*)'STABITY CRITERION ALERT: REDUCE TIME STEP'
END IF

! Total simulation time is 24 hours. How many interation steps is this?

ntot = 15.0*3600./dt  ! Remember that an hour has 60 minutes of 60 seconds each

! Data output at every hour of the simulation.

nout = 1.0*3600./dt

! Open file for output.
! The first (unit) number is a referene for output (see below).
! The "file" statement specifies the desired file name. 
! The statement "formatted" means ascii output.
! The status "unknown" implies new creation of a file if this does not exist,
! otherwise an existing file will be overwritten.
! Be careful not to overwrite files that you might need in the future.   


open(10, file = 'daimu_decay.txt', form = 'formatted', status = 'unknown')

! Write initial time and concentration to this file.

WRITE(10,*)0,100.0,100.0,100.0 

! Now comes the start of the iteration loop.
! It means that the statements between the DO-loop start and 
! the corresponding "END DO" statement are repeated for n = 1 to
! n = ntot at steps of 1. In this case, statements are repeated
! a total number of 864000 times. Try to do this on a piece of paper... 

!****** Start of iteration *******
DO n = 1,ntot   
!*********************************

CN_explicit = C_explicit*fac_explicit ! prediction for next time step
CN_implicit = C_implicit*fac_implicit ! prediction for next time step

time = REAL(n)*dt ! time counter

CTRUE = CZERO*exp(-kappa*time)  ! exact analytical solution

C_explicit = CN_explicit  ! updating for upcoming time step
C_implicit = CN_implicit  ! updating for upcoming time step
! Data output if the counter i is a constant integer multiple of nout.
! For instance, nout = 30 produces outputs at i = 30, 60. 90, ....
! Note that two equal signs (==) have to be used in the IF statement.

IF(mod(i,nout) == 0)THEN

   ! Output to unit 10. The associated file name is given above. The "*" creates an automated format.
   ! Here we produce three output columns: the first is the time in hours, the second our prediction, and the
   ! third the exact solution. IF-statements with more than 1 line have to use a "THEN" and an "END IF".

   WRITE(10,*)time/3600.0, C_explicit, C_implicit, CTRUE

   ! We can also write a note a message to the screen here.

   write(6,*)"Data output at time = ",time/3600.0," hours"

END IF

!****** End of iteration *******
END DO
!*******************************

STOP
END PROGRAM decay

import numpy as np
import matplotlib.pyplot as plt

plt.figure('Decay Problem')

data=np.loadtxt('E:\公众号\Ocean Modelling for Beginners\第二章\daimu_decay.txt').T;
res_explicit = data[1,:]
res_implicit = data[2,:]
analytical = data[3,:]

time = np.arange(0, 16, 1)
plt.plot(time, analytical, color='black', linestyle='dashed', label='analytical')
plt.plot(time, res_explicit, label='explicit')
plt.plot(time, res_implicit, label='implicit')

plt.legend()
plt.title('Decay Problem')
plt.xlabel('Time (hours)')
plt.ylabel('Concentration (100%)')
plt.ylim(0, 100)
plt.xlim(0, 15)
plt.xticks(time)
plt.grid(True)

plt.show()