猴子也能看懂的海洋数值模拟（一）——The Decay Problem

& I) y. S$ [% r/ V8 m( D4 |9 W

一、内容梗概

“Ocean Modelling for Beginners”第一章主要是一些编程语言和软件的介绍还有一些序言之类的内容，在此就不多做解释了。“The Decay Problem”是书中第二章最开始引入的一个问题，借此向读者介绍数值计算的基本知识。

{8 W7 J1 G5 H* M

“The Decay Problem” 的中文怎么称呼都好啦，衰退问题、腐败问题，总之就是假设某种物质的浓度是 C0C_0 ，这种物质的浓度随着时间变化会按照一定的速率逐渐减少，例如酒精挥发等等，这样一个过程可以由下面一个方程来表示：

dCdt=−κ⋅C\frac{dC}{dt} = - \kappa \cdot C \\

这是一个一阶常微分方程，其中变量 CC 表示某种物质的浓度， tt 代表时间， κ\kappa (希腊字母：kappa)是一个常量，影响衰减的速率。有了这样一个式子，我们就可以求解出任意时刻 tt 这种物质的浓度 C(t)C(t) 了。这个方程解起来也比较简单，对等式两边同时积分就能得到解，也就是未知函数 C(t)C(t) :

C(t)=c⋅e−κ⋅tC(t) = c \cdot e^{-\kappa \cdot t} \\

然而，由上式可以看出这个解不是唯一的(注意e前面的c指的是任意常数不是浓度)，要确定唯一的解，还要在一些自变量点上给出 C(t)C(t) 的值，也就是边界条件。这里我们可以设 t0t_0 时刻物质的浓度 C(t0)=C0C(t_{0}) = C_0 ：

' f# _. S, L5 }

{dCdt=−κ⋅C,t≥t0C(t0)=C0\left\{ \begin{array}{l} \frac{dC}{dt}= -\kappa \cdot C, t \geq t_0 \\ C(t_{0})=C_0 \end{array} \right. \\

在合理假定下，我们就可以认为从 t0t_0 时刻对应的初始浓度 C0C_0 开始，上式中的常微分方程决定了未知函数 C(t)C(t) 在 t_0">t>t0t > t_0 时的变化情况，也就是说这个边界条件可以确定常微分方程的唯一解。相应地，C(t0)=C0C(t_{0}) = C_0叫做初始条件，像上式一样带初始条件的常微分方程问题被称为初值问题。

那么讲了这么多，该如何利用计算机来求解这样的问题呢？这里就需要介绍另一个概念——有限差分(finite difference)了，形如 f(x+b)−f(x+a)f(x+b)-f(x+a) 的式子称为有限差分，将上式除以 b−ab-a ，我们就得到了差商，在数值计算中用有限差分来近似替代微分是一种非常重要的方法，有限差分通常分为三种前向差分(forward difference)、后向差分(backward difference)和中心差分(central difference)。前向差分的形式为 f(x+h)−f(x)f(x+h) - f(x) ，其中 hh 是步长，假设步长为 Δt\Delta t ，我们可以将上面的常微分方程改写成前向差分的形式：

dCdt=C(tn+Δt)−C(tn)Δt=−κ⋅C(tn)\frac{dC}{dt} = \frac{C(t_n + \Delta t) - C(t_n)}{\Delta t} = -\kappa \cdot C(t_n) \\

/ K1 G' N8 A" E ^$ D

换一种表示方式，让 C(tn)C(t_n) 为 CnC_n 也就是 tnt_n 时刻的浓度为 CnC_n ，则 tn+1t_{n+1} 时刻的物质浓度为C(tn+Δt)C(t_n + \Delta t) 即 Cn+1C_{n+1} ：

Cn+1−CnΔt=−κCn\frac{C_{n+1} - C_n}{\Delta t} = -\kappa C_n \\

( x+ o# B1 \" q! G( l

整理一下左右两边的项可以得到：

( ~6 f0 Z* O& M. @% L' h# ?' N

Cn+1=Cn−Δt⋅κ⋅Cn=(1−Δt⋅κ)CnC_{n+1} = C_n - \Delta t \cdot \kappa \cdot C_n = (1-\Delta t \cdot \kappa) C_n \\

这样的有限差分格式叫做显式格式(explicit scheme)，利用上式，我们就可以从t0t_0开始“步进式”的计算出后续时刻的函数近似值，在计算机程序中利用一个循环就可以完成。细心的同学可能注意到在这种方法下 Δt\Delta t 和 κ\kappa 的值是不能随便取的，试想一下如果 1">Δt⋅κ>1\Delta t \cdot \kappa > 1 那么 t0t_0 时刻之后的浓度就会变成负数，现实中根本不可能存在这种情况。因此我们的时间步长 Δt\Delta t 和 κ\kappa 需要满足下面一个条件：

) G) ?. Y/ ]6 |1 W/ \

% {) D! m- B2 f4 ^

该条件被称为数值稳定性条件(condition of numerical stability)。换句话说，只有当上式被满足时，我们的预测值，也就是通过有限差分发计算出的t时刻的物质浓度，是稳定的。可以看出这种方法下， Δt\Delta t 的最大值要受 κ\kappa 的限制。

除了上述的显式格式以外，还有两种格式，分别是隐式格式(implicit scheme)和混合格式(hybrid scheme)。

, m5 g# E( A: |

隐式格式

这种格式是通过后向差分得来的，后向差分其实就是用 xx 和 x−hx - h 替代显式差分中的 x+hx + h 和 xx ，将之前的微分方程写成后向差分的形式：

Cn+1−CnΔt=−κCn+1\frac{C_{n+1} - C_n}{\Delta t} = -\kappa C_{n+1} \\

0 G* `/ o5 K( D. z. x7 W; |2 I

和之前一样整理下等式左右两边，得到了这样一个结果：

$ M% N" E! i$ G: ]) p) A4 l

Cn+1=Cn(1+Δt⋅κ)C_{n+1} = \frac{C_n}{(1 + \Delta t \cdot \kappa)} \\

与显式格式相比，一个很大的好处是 Δt\Delta t 的取值不再受到 κ\kappa 的限制了。而且分母始终是大于1的( 0" >κ>0\kappa >0 )，因此浓度一定会随着时间逐渐减少。

4 u6 \9 y% D: C- \4 E, m

混合格式

如字面意思，这种方法混合了上述两种格式，形式如下：

2 ~! g9 j6 V# T9 Q

Cn+1−CnΔt=−α⋅κ⋅Cn+1−(1−α)⋅κ⋅Cn\frac{C_{n+1} - C_n}{\Delta t} = - \alpha \cdot \kappa \cdot C_{n+1} - (1- \alpha) \cdot \kappa \cdot C_n \\

整理一下可以得到：

# b4 P+ d1 M2 a1 q: d/ v+ J

Cn+1=1−(1−α)⋅κ⋅Δt1+α⋅κ⋅ΔtCnC_{n+1} = \frac{1-(1-\alpha) \cdot \kappa \cdot \Delta t}{ 1 + \alpha \cdot \kappa \cdot \Delta t} C_n \\

其中 α\alpha (希腊字母：alpha)是一个权值，决定哪一种格式所占的比重大，当 alpha=1alpha =1 时这种方法完全变成了隐式格式，而当 α=0.5\alpha = 0.5 时，这种格式被称做半隐式格式(semi-implicit scheme)。

/ U/ t' j0 e4 z, K! w

二、参考程序

9 d2 ~2 J7 T. t6 z! X$ ]

上面就是“The Decay Problem”练习之前的大致内容和基本概念，接下来就进入编码阶段了。题目的要求是： κ=0.0001s−1\kappa = 0.0001 \quad s^{-1} ， Δt=3600s\Delta t = 3600 s ，初始浓度 C0C_0 为100(%)，模拟15个小时的浓度变化。 程序中我将三种方法都写了出来方便后面比较它们的异同，其中混合格式的 α\alpha 值取了0.5。同时程序中还包含了一个求 tt 时刻真实值的方法，利用的是文章一开始积分求出的 C(t)C(t) ， c=C0c = C_0 。和猴子也能看懂的海洋数值模拟（零）里面说的一样，计算的部分用C++来写，数据处理和绘图利用python来完成，两部分代码如下：

e9 j" E' r2 C1 d% a0 [! s% ]

数值模拟程序

#include <fstream>#include <cmath>#include <vector>#include <functional> const double kappa = 0.0001f; const double alpha = 0.5f; ' v/ b% o# I, Q: E0 p2 \4 i! I" x- I const double dt = 3600; // 3600 seconds // lambda expression auto explicit_scheme = [](double conc)->double{ return (1 - dt * kappa) * conc; }; auto implicit_scheme = [](double conc)->double{ return conc / (1 + dt * kappa); }; & z% P$ V) `6 N5 t2 V auto hybrid_scheme = [](double conc)->double{ return (1 - (1 - alpha) * kappa * dt) * conc / (1 + alpha * kappa * dt); }; ' F* {; S3 ^8 J- b void Decay( double concentration, 9 p1 k1 J1 c7 x8 o7 L9 [2 Z int steps, const std::function<double(double)>& scheme, # P# K6 r( {- l1 f+ h3 m4 q8 N const char* result_path ) ! `/ _2 q0 F. Y; l4 i. L; [ { 9 C; i+ V; v4 V7 r i. ^4 H8 i std::vector<double> res{concentration}; for(int i = 0; i < steps; i++) # f8 c7 `% g% c) Z& [7 o$ l* Z6 }- D { # I3 Q% r% O+ j double conc1 = scheme(res[i]); res.push_back(conc1); - F9 s2 u6 T/ ^5 z6 a7 t! `4 e } ! }3 X2 J! `4 Z9 C' n: S # v! L% \. p. G% j5 p std::ofstream file(result_path); if(file.is_open()) ! q( d# R. j$ p% b { + e7 M; ]( [5 {8 i for(auto val : res) 3 G. B2 Z/ O8 ^; C# q6 c. x0 M file << val << "\n"; file.close(); % a0 t9 @5 c8 B6 U } 1 U- h% \; ?2 m3 I1 w; G } . Q" v# @+ v3 q! |9 W 5 \: S# Z& t$ I, w, m0 E8 R5 K void AnalyticalSolution(double concentration, int steps) 6 V9 B' O$ B0 g3 N { std::vector<double> res{concentration}; ' I! x1 j! |$ y4 I" z3 v double t = 0; for(int i = 1; i <= steps; i++) % l: M, T% R) {# m { t += dt; double ct = std::exp(-kappa * t) * concentration; res.push_back(ct); } ; N2 j6 [! \: t! Q1 z8 F std::ofstream file("res_analytical.txt"); if(file.is_open()) { 8 ]3 t' V/ ]( R for(auto val : res) file << val << "\n"; file.close(); } } int main() { / t% y9 X' ^% _9 E // the initial concentration is 100 % ( b1 A6 P* B3 S double concentration = 100; / F9 O$ J: g0 t* A // which means the total simulating time is steps * dt seconds ) C0 F" @; e5 U* g/ m* T int steps = 15; Decay(concentration, steps, explicit_scheme, "res_explicit.txt"); Decay(concentration, steps, implicit_scheme, "res_implicit.txt"); ' }7 `& Z; n$ F' E) T: T7 f: ^% j( c Decay(concentration, steps, hybrid_scheme, "res_hybrid.txt"); AnalyticalSolution(concentration, steps); 4 O$ m1 [- k" i* w return 0; " A4 e3 t+ |9 I. _ }

绘图程序

import numpy as np 7 y& r6 a- o! |( l: ] ?8 z q4 t import matplotlib.pyplot as plt + d/ M" o5 M3 \, b1 O) V plt.figure(Decay Problem) 6 b1 s* M4 }% V res_explicit = np.loadtxt(res_explicit.txt).T res_implicit = np.loadtxt(res_implicit.txt).T res_hybrid = np.loadtxt(res_hybrid.txt).T analytical = np.loadtxt(res_analytical.txt).T + r k7 m, o6 h' D+ I7 J& L1 \ time = np.arange(0, 16, 1) plt.plot(time, analytical, color=black, linestyle=dashed, label=analytical) + l; Y9 [8 i" }1 B+ d; e N* v( W plt.plot(time, res_explicit, label=explicit) 2 g4 s, e# g8 x0 u, Q0 @ a( Q; c( w% j plt.plot(time, res_implicit, label=implicit) plt.plot(time, res_hybrid, label=hybrid) ; L9 X! e0 s; [ plt.legend() plt.title(Decay Problem) plt.xlabel(Time (hours)) plt.ylabel(Concentration (100%)) : ]* H: L2 [/ O( \% B. g O) t plt.ylim(0, 100) plt.xlim(0, 15) plt.xticks(time) & u- I0 [0 R0 C; d plt.grid(True) % q+ \0 _5 g6 R, V' V plt.show()

三、运行结果

结果如下图1所示，程序模拟了初始浓度为100%的某物质从零时刻开始15个小时的浓度变化情况，分别采用了显式、隐式以及半隐式的差分方法，对比黑色虚线的真实值可以看出同一时刻显示差分法的计算出的浓度比真实值要少，相反，隐式差分法要多一些。令人惊讶的是半隐式差分法（hybrid），几乎和真实值的曲线吻合，从此图可以直观的看出三种不同的差分方法的区别。

3 g" K3 A2 o4 }# t

图 1. Decay Problem程序模拟结果

参考文献

Jochen Kampf. Ocean Modelling for Beginners. Springer. 2009. 喻文健. 数值分析与算法(第2版). 清华大学出版社. 2015.Wikipedia. 有限差分. https://en.wikipedia.org/wiki/Finite_difference.

之前的文章

猴子也能看懂的海洋数值模拟（零）

+ `5 O8 R; X( @8 I+ f $ K4 u# T6 x, b3 o