《计算物理》PPT课件

第八章常微分方程的数值解法,工程技术和自然科学中的许多实际问题，在数学上往往可以归结为求解微分线性方程的形式。高数中我们学过的解析方法主要是用在一些简单和特殊的微分方程求解中，而对于大量的一般形式的微分方程，不能用解析方法求出其精确解，而只能用近似方法求近似解。,在具体求解微分方程时，要有某种定解条件，微分方程与定解条件在一起称为定解问题。定解条件有两种：一种是给出积分曲线在初始点的状态，称作初始条件，相应的定解问题称作初值问题；另一种是给出积分曲线首尾两端的状态，称作边界条件，相应的定解问题称作边值问题。,本章主要研究的是一阶常微分方程初值问题的数值解法。,一阶常微分方程的初值问题是指形如,的问题。当函数f(x，y)满足,时初值问题（8-1）的解y=y(x)存在且唯一，此条件称为李普希兹(Lipschitz)条件。L为正常数，称为Lipschitz常数。,（8-1）,我们讨论的数值解法基本思想如下：在区间a，b上插入一系列离散点,寻求未知函数y(x)在这一系列离散点上的近似值yi(i=1,2,n)，而yi称为初值问题的数值解。常微分方程的数值解法期待建立y(xi)的近似值yi的递推格式，通过适当选取的初始值，求得y(x)在各点上的近似值。,第二节欧拉方法,欧拉法是解常微分方程组初值问题最常用的方法，我们用它来说明求解微分方程的基本的技巧和概念。,欧拉法的计算格式是,(8-2),这是一个递推公式，用它可以在已知yi时推出yi+1，下面给出公式的几种分析解释：,（1）将微分方程中出现的导数,用差商,近似代替，并以yi和yi+1分别,表示y(xi)和y(xi+1)的近似值，就得到,解出yi+1即可得到（8-2）式。h=xi-xi-1称为步长；一般取等步长h=(b-a)/n.,（2）从xi到xi+1对式,积分，得,（8-3）,如果采用左矩形公式近似计算上式中的积分，并以yi和yi+1分别表示y(xi)和y(xi+1)的近似值，也可得到,称为向前欧拉公式，也称为显式欧拉公式。此公式对给定的y0，可以算出y1，然后一步一步算出yi+1。,（8-2）,如果对（8-2）式中的积分采用右矩形公式，将得到计算格式,称为向后欧拉公式，也称为隐式欧拉公式。由于未知量yi+1含在右端，因而不能直接计算，一般使用迭代法进行求解。,(8-4),（3）利用y(x)在x=xi处的一阶泰勒公式，得,舍去余项，并以yi和yi+1分别表示y(xi)和y(xi+1)的近似值，就得到,梯形公式,上述的两种公式是用矩形求积公式来计算（8-3）式右端的积分而得到的，由于矩形公式精度较低，所以导出的欧拉公式精度也很低。如果用数值积分的梯形公式计算（8-3）式右端的积分，则有,（8-5）,这就是梯形公式。梯形公式由于其右端含有未知的yi+1，所以也是隐式方程。实际中也使用迭代法进行求解。,以f(x0，y0)为斜率，通过点(x0，y0)做一条直线，它与直线x=x1的交点就是y1。依此类推，yi+1是以f(xi，yi)为斜率过点（xi，yi）的直线与直线x=xi+1的交点，所以欧拉法也称为欧拉折线法。,欧拉方法的几何意义,预估校正法(改进欧拉法),显式欧拉公式计算工作量小，但精度低；梯形公式虽然提高了精度，但为隐式并且每一步都要解一个方程，计算量较大。基于这种想法，在实际计算中我们将这二种公式结合起来使用。,具体做法是：先用欧拉公式（8-2）求出一个初步的近似值，记作,，称之为预报值；,然后使用梯形公式作校正，即用预报值,代,替（8-5）右端的yi+1，再直接计算，得到的值yi+1就是校正值，这样建立起来的预估校正公式称作为改进欧拉公式，即,预报,校正,写成迭代形式,或表示为（编程时使用）,欧拉公式的截断误差与精度,局部截断误差：对于数值方法，若假定在前n个节点出的数值解等于其精确解，即yi=y(xi)(n=0,1,2，,n-1),记Ri+1=yi+1-y(xi+1),称Ri+1为数值方法在节点xi+1处的局部截断误差。局部截断误差在一定程度上反映了该方法的精度。通常用泰勒展开式来讨论数值方法的局部截断误差。,下面讨论欧拉法（以向前欧拉公式为例）的局部截断误差，将泰勒展开式,与欧拉公式,对照可得,定义若一种数值方法的局部截断误差为,，则称该方法具有P阶精度或,P阶方法。,由推导得欧拉公式具有一阶精度，故称一阶方法；而梯形公式有二阶精度，称二阶方法。,例1用向前欧拉公式和预估校正公式求解初值问题并比较它们的精度,取h=0.1；方程的准确解是,解：（1）用欧拉法计算,显然,x0=a=0，n=10，b=1，y0=1,由向前欧拉公式,（2）用预估校正公式计算,用以上公式计算公式可以求出各结点上的数值解，计算结果见下表。,同精确解比较，第2列欧拉公式的结果大约只有两位有效数字，而用预估校正公式的结果则有三位有效数字。,例2用欧拉法和预估校正公式解初值问题：,解：取f(x)=y2，y0=1，h=0.1,（1）用欧拉法求解,该方程的精确解是,，计算结果如下表。,（2）用预估校正公式求解用下面的迭代公式，对每个迭代4次，k=1，2，3，4。,计算结果如下表。,两种方法进行比较，发现用预估校正公式误差较小。,第二节龙格库塔法,在上一节中，我们得到了一些求微分方程近似解的数值方法，这些方法的局部截断误差较大，精度较低，我们希望得到有更高阶精度的方法。,考察差商，由微分中值定理，存在点，使得,其中称为y(x)在xi，xi+1上的平均斜率。这样，只要对k*提供一种近似算法，便可以导出一种计算格式，这就是龙格库塔方法的基本思想。,（8-6）,即:,一阶龙格库塔方法,如果以y(x)在xn处的斜率作为y(x)在xi，xi+1上的平均斜率k*即,则得,k*,这就是欧拉法,二阶龙格库塔方法,考虑在xi，xi+1上取两点xi，xi+p（0p1）的斜率值k1，k2的线性组合1k1+2k2作为k*的近似值(1、2为待定常数)，此公式一般形式可写成,其中k1=f(xi，yi)，k2为xi，xi+1内任意一点xi+p=xi+ph(0p1)的斜率f(xi+p，y(xi+p)。,由于y(xi+p)并没有给出，仿照改进欧拉公式的构造思想，得到下面的公式,（8-7）,这样构造出的公式为,k1,k2,k*,公式中含有三个参数1，2和p，如果我们适当选取参数的值，可以使公式的局部截断误差为O(h3)。,对k1和k2作泰勒展开,代入（8-7）得,（*）,又y(x)在xi处的二阶泰勒展开式为,当x=xi+1时，,，有,（*）,（*）,比较（*）与（*）的系数即可发现，要使公式（8-7）的局部截断误差满足,，即要求公式具有二阶,精度只要下列条件成立即可。,（8-8）,满足条件（8-8）的一簇公式统称为二阶龙格库塔公式。,特别的，当,塔公式就成为改进欧拉公式。,时，龙格-库,改进欧拉公式就是以y(x)在xi和xi+1处的斜率k1和k2的算术平均值作为y(x)在xi，xi+1上的平均斜率k*来进行计算的。,如果取,龙格-库塔公式就称为变形的欧拉公式，其形式为,（8-9）,此处的,就是欧拉方法预报出的中点,处的近似解；而,等于中点的斜率值,则近似,；,时，,所以公式可以看作用中点斜率近似代替平均斜率k*，因此，公式（8-9）也称作中点公式。,三阶龙格库塔公式,为了提高精度，考虑在xi，xi+1上取三点xi，xi+p，xi+q的斜率值分别为k1，k2，k3，将k1，k2，k3的线性组合作为平均斜率k*的近似值，其中,xi+q=xi+qh（0pq1）,这时的计算公式为,其中,为了预报点xi+q的斜率k3，一种很自然的想法是用欧拉法预报，即取,但是，这样做效率比较差。因为在区间xi，xi+q内已有两个斜率可以使用，所以把k1，k2的线性组合作为xi，xi+q上平均斜率的近似值，当然比用欧拉法预报y(xi+q)精度要好。由此，得到y(xi+q)的预报值。,于是可取,从而得到三点计算公式的形式,（8-10）,下面列出其中的一种形式,类似前面的二阶龙格-库塔公式的推导，利用泰勒展开法适当选择参数1，2和p，q，r，s的值可以使上述公式具有三阶精度，即局部截断误差为O(h4)。这一类公式统称为三阶龙格-库塔公式。,（8-11）,高阶龙格-库塔法,为了进一步提高精度，在xi，xi+1上可取多个点，预报相应点的斜率值，对这些斜率值加权平均作为平均斜率值。利用泰勒展开，比较相应系数，从而确定在尽可能高的精度下有关参数应满足的条件。,从理论上讲，可以构造任意高阶的龙格-库塔公式，但实践证明，高于四阶的龙格-库塔公式，不但计算量大，而且精确度并不一定提高。在实际计算中，四阶龙格-库塔公式是精度及计算量较理想的公式。,经典的四阶龙格-库塔公式为,（8-12）,龙格-库塔法的步长选择,上面所讨论的龙格-库塔方法是定步长的，单从每一步来说，步长h越小，局部截断误差就越小；但随着步长的缩小，不但引起计算量的增加，而且也有可能引起舍入误差的严重积累；而步长h太大又不能达到预期的精度要求，所以怎样选择合适的步长h，这在实际计算中是很重要的。,我们以四阶龙格-库塔公式为例，从节点xi出发，先以某个h为步长求出一个近似值，记为yi+1(h)，然后将步长折半，即取h/2为步长从xi计算两步到xi+1，再求得一个近似值yi+1(h/2)。,通过计算|yi+1(h/2)-yi+1(h)|是否成立，成立则步长选择h/2，否则继续将区间分成h/4，再次重复上述步骤至满足精度要求。,例使用高阶R-K方法计算初值问题,解：,(1)使用三阶R-K方法,其余结果如下:,(2)如果使用四阶R-K方法,ixik1k2k3yi1.00000.10001.00001.10251.25551.11112.00000.20001.23451.37551.59451.24993.00000.30001.56241.76372.09221.42844.00000.40002.04042.34232.86581.66645.00000.50002.77683.25874.16341.9993,其余结果如下:,ixik1k2k3k4yi1.00000.10001.0