数学物理方程爆破解的数值诊断方法

第1章一阶常微分方程柯西问题爆破解的诊断 在这一章，我们通过一个具体的例子来介绍我们所关心的数学问题.考虑如下一阶常微分方程的柯西问题(它是化学动力学中一个简单的热传导模型): (1.0.1) 为了确定起见，我们代入以下参数: (1.0.2) 可以验证，函数 (1.0.3) 是问题(1.0.1)–(1.0.2)的一个解.我们将函数(1.0.3)(其中参数u0的取值见(1.0.2))在图1.1中画出.可以看出，函数(1.0.3)在有限的时间内(在t=1内)趋向了无穷，并且在t=1点没有定义.但是我们认定柯西问题的解是定义在连通集合(一段区间)内的函数，所以问题(1.0.1)的解是函数(1.0.3)限制于t∈[0，1)时的部分.对于此类情形，我们说方程的解经历了爆破(blow-up). 需要注意的是，在大部分自然科学领域的实际问题中，某个参数在有限的时间里趋向正无穷的情形是不可能存在的.而这个参数可能在相应的简化的数学模型中的解中趋向无穷.在这种情况下，我们可以认为真实问题超出了简化数学模型的可应用范围.例如，在化学反应的热传导问题中，只有在无限多的反应物集中在相当小的体积的假设下，才有可能在t=1附近达到无穷.显然，这种情况在实际中是不会遇到的，因为反应物的量总是有限的，这会*终导致反应的停止，以及在某个固定的小于1的时间点t时停止放热(在这种情况下，总放热量是有限的).因此，如果我们弄清所考虑物理问题的合理模型，并同时考虑到守恒定律和反应传播速率，那么所得的数学模型将给出从起始时间点开始都存在的解.但是，在实际情形中精确模型的构造(由于物理过程的复杂性，缺乏准确的数据等等)或者甚至是数值解的求解(由于方程式的极端复杂性)都可能是无法实现的.因此，它们通常仅限于描述过程特征的简化模型.近似模型中的解趋向于无穷可能对应真实的爆炸(非常大的热量释放). 但即使简化后的模型也十分复杂，以至于解决它们的唯一有效方法就是使用数值方法获得近似数值解.对于任意t的取值都存在数值解(对于适当近似值的选择可以查阅1.2节的备注2).我们以研究过的标准模型(1.0.1)–(1.0.2)为例进行演示. 图1.1满足方程(1.0.1)和参数集合(1.0.2)的函数. 在阴影线表示的区域里函数u(t)不是相应的柯西问题的解 1.1数值解的搜寻 将问题(1.0.1)改写为以下形式 (1.1.1) 此处f(u)=u2. 引入对于时间t的等长区间TM，间隔大小为.注意到，这个区间有M+1个节点，或者等效地，M个小间隔.在网格TM的节点处引入函数u(t)的网格值. 对于自治系统(1.1.1)写出单阶段Rosenbrock算法群ROS1[16，22–24] 这里w1定义于方程 (1.1.2) 这里a11是定义算法性质的参数(关于这一点的详细内容请看文献[16，Chapter1，Section1.2.3]).特别地，对于a11=0，算法(1.1.2)退化为一个明确的单阶段Runge-Kutta算法ERK1(欧拉算法);当a11=1时，退化为反欧拉算法DIRK1;当a11=1/2时，退化为“半隐式”算法;当a11=(1+i)/2时，退化为带复系数的单阶段Rosenbrock算法CROS1. 以下是MatLab函数的示例，该函数根据算法(1.1.2)实现对问题(1.1.1)的数值求解. 说明需要注意的是，当访问向量u和t的分量时，所有的索引需要位移+1位(与上述解析公式相比)，因为在MatLab中数组元素的编号从1开始(因此). 下面，我们给出可以运行函数“ODESolving”和展示结果的代码. 上述程序对参数a11进行不同的取值，其结果(对测试问题(1.0.1)–(1.0.2)进行求解)展示在图1.2中. 由图1.2可以发现，对于所有的取值都存在一个数值解.因此出现了一个自然的问题:如何分析确定解爆破的事实并确定哪个数值解是可信的，而哪些不能？对于那些无法求出解析解的复杂方程，这个问题尤为重要.接下来我们将详细研究这个问题. 图1.2在相同的时间网格M=50运行算法(1.1.2)，以及通过对参数a11进行不同的取值，对测试问题(1.0.1)–(1.0.2)求解的结果 1.2解的爆破现象的数值分析 我们在此书中讨论的解的爆破现象的数值分析方法基于误差的后验渐近精确估计的计算(见[16，第2章]). 假设我们通过使用对时间t均分网格TM，步长为的p阶精确算法已经找到了柯西问题(1.0.1)的网格数值解.这意味着，对于所有的节点t∈TM都成立等式 (1.2.1) 这里的上标(M)表示对应于M个时间网格的数值解u(M)(t). 现在我们将等式(1.2.1)写成以下形式 (1.2.2) 在上式中，我们将近似解u(t)的误差u(t)-u(M)(t)的泰勒级数的首项提取出来.这里我们假设存在相应的连续导数. 假设我们通过已知p阶精度的算法得到u(M)(t)，t∈TM(注意到区间的个数M唯一确定了网格步长τ).这样，方程(1.2.2)含有两个未知量:u(t)和c(t).为了求出这两个未知量，我们需要两个方程.在r倍紧密的网格(就是说在含有rM个区间的网格上进行计算)，我们可以得到第二个方程: (1.2.3) 从方程(1.2.3)中减去方程(1.2.2)可以得到c(t)的表达式: 因此， (1.2.4) 等式(1.2.4)右端的分式被称为精确解u(t)泰勒级数的p阶项的后验渐近估计.它也是近似解u(t)的误差u(t)-u(rM)(t)的泰勒级数的首项的后验渐近估计. 因此， 从这里可以得出，当τ→0时，R(rM)(t)大于误差量u(t)-u(rM)(t)的泰勒级数中的所有其他项之和.这一项可以看作u(rM)(t)的误差的后验渐近估计.为简明起见，本书将忽略O(τp+1)(但此时不要忘记相应公式的渐近特征).即，此后我们将R(rM)(t)当作 (1.2.5) 这就是经典的龙格–龙贝格(Runge-Romberg)公式. 我们再在一个新网格上进行计算(网格有r2M个区间)，可以得到 注意到 (1.2.6) 通过这两个式子的范数的关系求出精确度有效阶数的表达式 (1.2.7)