最近网络疯传疫情模拟预测结果。CFD界在看到这个模型之后,感觉很欣慰。数值模拟手段正在走进千家万户。
另一方面,网络上流传的各种模型,存在一些缺陷。其实类似的问题在计算流体力学(CFD)中特别常见,比如大名鼎鼎的群体平衡模型Population Balance Equation(PBE)
Population是什么,是人口。Population Balance就是人口平衡。
PBE在CFD这面求解已经几十年了。只不过从来没和医学交叉过。
在这个小文中,CFD界简单讨论,如何从偏微分方程的角度,更完善的描述感染人口的扩散。
更完善的模型,需要和医学领域,交叉合作研究。
我们认为P表示感染人口数量,其为一个随时间、空间变化的函数,因此,做导数的时候,为偏导!:
对P对时间做微分,表示感染人口数量的变化。最简单的情况:
如果S=0,则P不增不减,表示新增感染人口与治愈人口数量相等。这个就是最简单的网络上疯传的SIP模型,可以看出,这个是算法求解中最简单的数值模型。
这只是一个ODE系统,由于疫情可控,方程不存在任何刚性,小学生都能求解。
1. 医学与流体力学交叉后的模型
对于我们CFDer,搞偏微分方程的,咱们玩点复杂的。
我们考虑感染人口的传播特性,在这里考虑传输人口的扩散行为,与主动对流行为,扩散行为可以表示为:
在这里Γ表示感染人口的扩散性传输。Γ为一个二阶张量!感染人口的主动对流行为可以表示为:
其中的U表示感染人口的主动传播方向,为一个矢量。将这些项放在一起,构成
感染人口传输方程:
方程未知量是P(感染人口),为随着空间,时间变化的量。
求解这个方程,可以得出感染人口的时间、空间分布。也就是这张疫情地图是可以计算机模拟出来了的。
在这里暂且将它称之为感染人口传输方程。这个方程强烈区别于网络流传的SIP模型以及前几天刷屏的B站模型。
在这里要强调,在这里没有调用任何的模型去封闭。因此感染人口传输方程是完全精准的。
感染人口传输方程还不能求解,因为一个方程存在多个未知量。
这个方程中,除了P是未知变量,U、Γ、S1、S2均为未知变量。
下面我们对方程进行封闭。
2. 模型封闭
首先是U,U表示感染人口的主动传输,也就是这些感染人口,倾向于去哪。在这里要对其进行模型化处理。
目前正值春运,我们认为U应该具备大城市间的主动传播特性,比如武汉- 北京方向的U的数值应该比较大,相反,全国各地到武汉的速度数值应该比较小。
类似的,经过这种模型化处理,有了感染人口主动传输速度U。
下面要对Γ进行封闭。最简单的模型为将其处理为一个标量,即感染人口无方向的向四周传播。
如果考虑感染人口的扩散反向,Γ为二阶张量,也即感染人口最后的扩散可能会呈现一个椭圆的形状或者其他。
这样,感染人口传输方程左侧则封闭。
感染人口传输方程右侧的S1和S2分别表示源项和沉降项,这俩项也和空间、时间变化有关。
考虑正的S1源项,其表示感染人口数量增加的时间变化,对S1进行模型化,需要考虑:
武汉地区的S1数值较大
随着时间变化,S1应该增加
其他
考虑负的S2源项,其表示感染人口数量减少的时间变化,这一项,大体反应了治愈率。具备以下特征:
在治愈率较低的情况下,S2的数值非常小
大城市医疗设备完善,S2是一个比较大的负值
疫苗的产生、特效药的生产,都会导致S2绝对值的增加
其他
这样,
感染人口传输方程封闭,即可求解,获得随时间、空间位置变化的感染人口数量。
在这里需要强调,如果U和Γ和P有关,则构成非线性PDE系统,
1)需要线性化处理;
2)需要迭代计算;
这就是我们CFDer专门干的事。
一个精准的U、Γ、S1以及S2的模型化处理,会影响最终预测的感染人口数量。
好了,在这里CFD界,只是从流体力学的角度,简单讨论了一下针对武汉疫情传播的数值模型。
其实大约2星期前,CFD界就考虑过写这个模型出来,但总感觉大家对这个不感兴趣,全是方程,乱七八糟的也看不懂。
不过,一个学科如果能用数学来表示,那才是最严谨的。
以本次疫情为例,这是一次完美的学科交叉。
感染人口传输方程的具体求解,感兴趣的,算一下吧。
如果有医学领域的关注者,可以研究些如何对方程各项进行封闭,具体求解,可以留给CFDer们,比如下图就是CFDer典型的计算过程:
先手算!再植入计算机求解!
本文是CFD界凌晨5点在梦中琢磨的,感谢你认同这个模型,如果不认同,就当我凌晨做梦好了。
有医学口的科研单位感兴趣搞起来不?联合申请个基金,我只需要5块买点草稿纸。