第四章-常微分方程数值解 1

第四章常微分方程数值解,初值问题数值解,边值问题数值解,1常微分方程初值问题数值解,一、问题的提出一阶常微分方程初值问题的一般提法是，求函数，使其满足微分方程：,(4.1),和满足初值条件：,(4.2),对式(4.1)进行积分，利用式(4.2)的初始条件，可得,(4.3),如果函数的形式比较简单，适合于直接积分，由式(4.3)即可求得随变化的连续函数。这就是分析解，也称精确解。如果不能直接积分，则可通过对式(4.3)的数值积分求值。也就是一阶常微分方程初值问题的数值解。,为了进行数值解，需将时间区域离散化。取时间步长，它将连续的时间区域划分为一系列离散的时刻，,(4.4),则(4.3)可以写为,若与都已知，只要反复应用式(4.4)即可算得任意时刻的值。,公式(4.4)也可用图4.1中曲线来说明。,-由图4.1(a)中A点表示,-由图4.1(b)中曲线表示,-积分为曲线所对应的面积,-若能精确计算积分，就可算得，即为图4.1(a)中B点。,(a),图4.1,(b),图4.1Eular法和改进Eular法几何意义。,常微分方程(4.1)的数值解，就是围绕着对式(4.4)如何进行数值积分展开的。目前比较常用的方法为RungeKutta法。,二、Euler法Euler法是实现式(4.4)数值积分最简单的方法，计算公式：,(4.5),已知与，而后重复运用式(4.5)，即可求得。,由求表示在图上为：过A点作曲线的切线，其斜率，此切线与的垂线相交于B1点。它即代表的值。,比较式(4.4)与(4.5)，对照图4.1，不难看到，Euler法的实质是：用矩形面积（高，底）代替曲线积分所对应的面积。显然这会引进相当的误差。为了减小这种误差，又提出了改进Euler法。改进Euler法的基本思想是，用斜线对应的梯形面积来代替曲线的积分。即,相应于式(4.5)的计算公式改写为：,(4.6),(n=0,1,2),(4.7),这就是改进Euler法计算公式。由于用代替了，不再落在点上，而在它的附近。,用式(4.6)代替(4.5)进行计算当然会准确些，但等式右边含有这个未知量，它是无法求得的。这使式(4.6)实际上不能用。为此，将上式右边的改为，即用近似值代替精确值，式(4.6)改为：,式(4.7)和式(4.5)的任务，都是在已知的条件下，求的值。但改进Euler法与Euler法在计算上有一个区别。Euler法的计算公式(4.5)中，待求的量只在等式左边出现，等式右边的量都是已知的，这类情况通常称显示方式。,而计算公式(4.7)中除了在等式左边出现外，还在等式右边函数内出现，这样在已知的情况下，必须通过解方程才能获得。因此，在改进Euler法中就出现了在Euler法中未曾出现的新问题，即求解有关的隐式方程。如果这种隐式比较复杂，则是否可以解以及如何求解，都是尚待解决的问题。,计算方法原则上讲，可以通过迭代的方法求值，例如用Euler法求得初值:,而后代入式(4.7)进行迭代计算，即：,(4.8),可以证明，当步长取得比较小，经过多次迭代后，将收敛于。这种改进Euler法在计算精度上虽然比Euler法优越，但计算量大，关键是用迭代法解方程不知是迭代多少次，在实际计算中，为了简化，只迭代一次算得。这时计算公式为：,(4.9),(4.10),通常称这类公式为预报校正公式，式(4.10)为预报公式，它预报第一个值，式(4.9)为校正公式，由此算出的校正值。将上述计算公式写成更简单的形式，其中,(4.11),(4.12),(4.13),K1的几何意义前面已经讲过，为的面积。,一般不落在图4.1(b)过的曲线上，因为时刻对应的只是的近似值。根据，而是一个单调下降的函数据，我们可预料，在曲线的下方。,注意：,K2的几何意义从图4.1(b)中可看到，是的面积。对应于，在求得K1后即可求得此值。,由于改进后的Euler法用斜直线对应的面积（梯形）代替曲线所对应的面积，它的精度仍受相当的限制，为了提高计算的精确度，又提出了Runge-Kutta法。,四、Runge-Kutta法问题仍然是已知（认为它是精确的）求的值，希望能算得更精确些（比前面介绍的方法），而计算工作量又不大。,受改进Euler法的启发，Runge-Kutta法的基本思想是，用在四个不同点上的数值的某种线性组合来近似代替积分,具体地说，Runge-Kutta法的计算公式如下：,其中,(4.14),(4.15),(4.16),(4.17),(4.18),五、误差分析无论是Euler法，改进的Euler法，还是Runge-Kutta法，它们的计算结果与精确解的结果总是有偏差，也称误差。方法不同，误差也就不同。因此，在采用何种方法进行数值计算时，对各种方法可能引起的误差应有一个定的分析。,这里所说的误差，是在假设，即取得准确值的前提下，采用数值求解方法求得的，与精确解之间的偏差。,精确解可根据Taylor级数展开公式写出：,余项,对于Euler法，由式(4.5)得,（因为,(4.20),(4.19),式(4.19)与式(4.20)相减，得,余项,(4.21),式(4.21)说明了，Euler法的计算公式实际上是Taylor级数截去余部（包括二次项）所得到的，这截去的部分也就是Euler法的误差，为“截断误差”，该截断误差的量为。即选用Euler法作数值计算时，每一步运算引进的误差为步长平方的数量级。,改进Euler法：,截断误差估计如下：,将按二元函数在与处展开Taylor级数：,(4.11),(4.12),(4.13),余项，,，,将，代入上式，得：,余项,再将代入上式，得：,余项,(4.22),将式(4.22)代入式(4.13)，得,余项,再将式(4.12)与式(4.23)代入式(4.11)，得,余项,(4.23),(4.24),式（4.19）与式（4.24）相减，得的差值为，这个差值是两个余项的差，它的数量级等同于余项的数量级。这说明改进Euler法的截断误差为，比Euler法提高了一个数量级，计算精度提高了。用上述同样的方法，可推导出Runge-Kutta法的截断误差。四阶Runge-Kutta法的截断误差为，计算的精度更高了。,讨论：从以上分析中可看到，截断误差的大小与所用的计算公式有关，同时，还与步长的大小有关。越小，截断误差也越小。为了得到较高的计算精度，通常取较小的步长。但同时，由于步长取得小，计算的量就增加了，与此相应的计算过程的舍入误差就会增加，这种误差有时可能起主要作用，因此，在总的误差分析中，要权衡舍入误差与截断误差两方面，选取合适的步长。,六、一阶常微分方程组初值问题的数值解在了解了Runge-Kutta法求解一阶常微分方程之后，不难把这个方法应用于一阶常微分方程组。由两个方程所组成的一阶常微分方程组初值问题表述如下：,把式(4.14)、(4.15)、(4.16)、(4.17)、(4.18)类似地用到方程(4.25)与(4.26)中，就可得到求解方程组的Runge-Kuntta法计算公式。,(4.28),(4.29),其中,(4.30),对于超过两个方程所组成的一阶常微分方程组，仍可用Runge-Kutta法。其计算公式同二个方程组成的常微分方程组是类同的。,2常微分方程边界值问题数值解,本节将以二阶线性微分方程为例来讨论两点边值问题，即考虑形式：,(4.2.1a),该方程中不显含。其边值条件可分为下面三类。,（4.2.1b）,（4.2.1c）,第一边值条件：,第二边值条件：,第三边值条件,(4.2.1d),解线性边值问题的差分方法,2.1差分方程的建立,今就线性微分方程的情况来讨论,（4.2.1）,用分点:,将区间a,b划分为n等分,称为结点。现在我们把在中求解的问题化为在这些结点上求的近似值的问题。为此，我们将(4.2.1)中的变量进行离散，具体做法如下：对于内部节点，将二阶导数用二阶中心差商来表示，便得到关系式,将它代入(4.2.1)便得到在结点上所满足的关系式,其中.,上式中略去最后一项，我们就可得到方程（4.2.1）的近似差分方程,（4.2.2）,这是含有n+1个未知数的线性代数方程组，方程的个数为n-1.,要使方程组（4.2.2）有唯一解，还需要由边值条件:补充两个方程。,对于第一边值条件，直接就得到了另二个方程,于是，就得到第一边值问题的差分方程组：,（4.2.3）,对于第二及第三边值条件，由于包含了导数，故边值条件也必须用差商来近似表示。因为我们不能利用在区间a,b以外的点，所以在引进（即）或（即）的近似式时，就不能利有中心差商。如果只要求误差是，可以用最简单的近似公式,（4.2.4）,(4.2.5),（4.2.6）,但要使误差达到，则需利用Newton等距插值公式，可以得到近似公式,把这些近似式代入边值条件中，再与方程组（4.2.2）联立，就可以得到对应的差分方程组。,对于第三边值问题，采用近似式(2.5),(2.6)，就得到差分方程组,（4.2.7）,这样我们便通过离散化的过程将微分方程的边值问题化为一个差分方程的边值问题了。,2.2解差分方程组的追赶法我们先考虑第一边值问题的差分解法。把（4.2.3）中的消去可得,(4.2.8a),(4.2.8b),故可用一般的主元素消去法来解上述方程组，可以利用方程组的特点将步骤加以简化。先从（4.2.8-1）中解出y1,（4.2.9-1）,方程组（4.2.8）的系数矩阵为三对角的。又因,将（4.2.9-1）代入（4.2.8-2）并解出y2,（4.2.9-2）,一般可设,（4.2.9-i）,为了求得及的递推公式，可将（2.9-i）代入（4.2.8-i+1），解出,因此有:,（4.2.10）为的递推关系式，解出之后由（4.2.9）就能求得，从（4.2.9）式可以看出，应满足关系式,（4.2.10）,小结：,已知，故应假设，这样便可逐个地按下标从小到大的次序求得和,又因，故由（4.2.9）式可以按下标从大到小的次序逐个求得。,求的过程下标由小到大，称为追的过程，求的过程正好相反，下标由大到小称为赶的过程，这就是解差分方程的追赶法。,上述计算过程实际上与主元素消去法是一致的。利用公式(4.2.10)计算和的过程相当于消元过程，利用（4.2.9）计算的过程相当于回代过程。,可以证明，追赶法有误差在传播中不增长的优点。追赶法在计算过程中是稳定的，故它是解差分方程组的有效方法。,例：考虑边值问题,取，则区间0，1被分成5等分，为其内点在内点上，我们以二阶中心差商,来代替二阶导数，得差分方程,它是一个线性方程组.可将它化成,（4.2.11）,的形式，其中（相当于q=0,r=0）,由于边界条件：，所以（4.2.11）式得,由（4.2.12）,令，得,然后可逐个求得,再由（4.2.11）与,逐个求得,第三边值问题的差分解法,方程和边界条件为：,可类似地建立差分方程。,所不同的是需要采取新的分点：,在内点上的差分方程为,在边界上均以差商来代替导数，于是从边值条件便又得两个方程:,这样便得到n+2个方程的线性方程组，经过整理后可以写成,（4.2.13）,（4.2.14）,（4.2.15）,与第一边值问题一样，我们把上述方程组用Gauss主元素消去法化成一组递推公式:,（4.2.16）,其中,（4.2.17）,所不同的是现在是未知的，需将(4.2.16)中的,然后就可由（4.2.16）逐个求出可以证明计算过程是稳定的。,计算方法：,先由（4.2.17）解出及;,然后再由(4.2.16)求出;,与(4.2.15)式联立起来解出;,3试射法,试射法（shooting）可用来解二阶或高阶的线性或非线性常微方程。这个方法的实质在于把边值问题化为初值问题来解。此时可以采用各种初值问题的单步法或多步法进行求解。下面将以第一边值问题与第三边值问题为例来介绍试射法是怎样用的。,今考虑二阶方程的第一边值问题,(3.1),试射法的思想就是设法确定的值m使满足初值问题的解也满足另外一个边值条件。也就是要从微分方程（3.1）的经过点,而具有不同斜率的积分曲线中，去寻找一条经过点的曲线。首先我们可以根据经验，或对方程作定性分析，或按照实际存在的运动规律，选取一个斜率m1，我们就用这个斜率进行试算，即解初值问题,这样便得到一个解，如果或（为允许误差），则即为所求的解。否则，可根据与的差距来适当地将m1修改（例如取）。这时，再解初值问题,于是又可得到另一个解。仿前，如果或满足不等式，则即所求的解，否则再对m值作适当修改。,事实上，假定为初值问题,的解，显然，在为确定的情况下，终值是初始斜率m的的函数：,于是，问题就变成求m使,这是一个代数方程（线性或非线性的），但由于的具体表达式往往是不明确的，所以解这方程一般是困难的。最简单的办法是由及用线性插值法来求出新的m值,(3.2),当然，用线性插值的依据是不足的，可是却往往用它来调节初值，因为它是方便易行。值得指出的是，如果有更适当的插值公式可利用的话，那末就有可能使尝试的次数有效地减少。,例：解边值问题,我们取第一次尝试解在（0，0）点的斜率为，通过计算得到，它与相差尚远，显然不符合要求，接下去我们取得到。此时由于仍不满足要求，所以就用线性插值公式,得到，相应的。这离又近了一些。我们再从重新开始按比例取，求得，它大于1。再由作线性插值，得，通过计算得。由于此值已比较接近1，故可再按比例取，算得。,再按比例取，算得。再由作线性插值求出。此时算得它与所要求的只相差0.001，计算可到此结束。,上述具体步骤