常微分方程数值解课件

科学研究和工程实践中，有很多实际问题的数学模型都是微分方程。利用微分方程理论，我们可以研究它们的一些性质，对实际问题进行分析。但是，只有极少数特殊的方程有解析解。对于绝大部分的微分方程是没有办法求出它的解析解的。,常微分方程作为微分方程的基本类型之一，在自然界与工程界有很广泛的应用。很多问题的数学模型都可以归结为常微分方程的定解问题。很多偏微分方程问题，也可以化为常微分方程问题来近似求解。,常微分方程数值解,本章讨论常微分方程的数值解法,考虑一阶常微分方程的初值问题,只要f(x,y)在a,bR1上连续，且关于y满足Lipschitz条件，即存在与x,y无关的常数L使对任意xa,b，和y1，y2R1都有|f(x,y1)f(x,y2)|L|y1y2|成立，则上述问题存在唯一解。,虽然有些微分方程初值问题的解是存在而且唯一的，但是，要求出其精确解yy(x)往往是非常困难的。为此，我们想法计算出解函数y(x)在一系列离散节点a=x0x1xn=b处的近似值yiy(xi)(i=1,2,n).这就是微分方程数值解。,节点间距hi=xi+1xi为步长，通常采用等距节点，即取hi=h(常数)。,求微分方程数值解的基本思想是：在这些节点上采用离散化方法，（通常用数值积分、微分、泰勒展开等）将初值问题化成关于离散变量的相应问题。把这个相应问题的解yn作为y(xn)的近似值。这样求得的yn就是上述初值问题在节点xn上的数值解。一般说来，不同的离散化导致不同的方法。,微分方程数值解：计算微分方程解函数y(x)在节点a=x0x1xn=b处的近似值yiy(xi)(i=1,2,n).,一、欧拉(Euler)法与改进欧拉法,1.欧拉法：,亦称为欧拉折线法,x1,称为后退的欧拉格式,x1,例1用欧拉法求初值问题在区间0,0.10上的数值解:,当h=0.02时在区间0,0.10上的数值解。,方程真解：,当h=0.02时在区间0,0.10上的数值解。,微分方程数值解的精度,欧拉法的局部截断误差：,=O(h2)=O(h1+1)欧拉法的精度为1阶。,一阶方程的初值问题与积分方程,是等价的，当x=x1时，,借助于数值积分，求y(x1)的值,若用矩形公式,2梯形公式,用梯形公式,则有,Euler公式,于是有递推格式：,称之为梯形公式。这是一个隐式的计算公式，欲求yn+1需要解方程。,Euler公式是求微分方程数值解的很好的方法，它算法简单，易于计算，但Euler方法有一个弱点就是误差较大不能保证精度要求。梯形公式的弱点在于需要解一个方程。,因此常采用将简单的Euler公式与梯形公式相结合的方法：即预测-校正法：,3改进欧拉法,改进的Euler公式,改进的Euler公式还有几种不同形式：,取步长h=0.1，分别用Euler方法和改进的Euler方法，求微分方程初值问题,在区间0，1.5上的数值解。,例题,Euler方法：y0=1,yn+1=yn+hf(xn,yn).,y1=y(0.1)=y0+hf(x0,y0)=1+0.1(120/1)=1.1,y2=y(0.2)=1.191818,y11=y(1.0)=1.784778,改进的Euler方法：y0=1,局部截断误差,二元泰勒公式：,设z=f(x,y)在点(x0，y0)的某一邻域内连续且直到有n+1阶连续偏导数，(x0h，y0k)为此邻域内任一点，则有：,欧拉法的截断误差：,改进欧拉法的截断误差：,改进的欧拉法的精度为2阶。,几何解释,建立高精度的单步递推格式。,单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。,考察改进的欧拉法，可以将其改写为：,二、龙格-库塔法,首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得,Step1:将K2在(xi,yi)点作Taylor展开,Step2:将K2代入第1式，得到,Step3:将yi+1与y(xi+1)在xi点的泰勒展开作比较,要求，则必须有：,这里有3个未知数，2个方程。,存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。取,注意到此式就是改进的欧拉法。,则有,类似的有三阶Runge-Kutta公式：,问题:为获得更高的精度，应该如何进一步推广？,其中i(i=1,m)，i(i=2,m)和ij(i=2,m;j=1,i1)均为待定系数，确定这些系数的步骤与前面相似。,最常用为四级4阶经典龙格-库塔法,4阶龙格库塔法,截断误差阶为O(h5)。,取步长h=0.1，用四阶龙格库塔法解初值问题,y=x2y(0x1)y(0)=1的数值解,解：初值问题的精确解为y=x2-2x+2-e-x,例题,用四阶龙格库塔公式,xiyiy(xi)(精确解)0.000000,y0=1.000000y(x0)=1.0000000.100000,y1=0.905163y(x1)=0.9051630.200000,y2=0.821269y(x2)=0.8212690.300000,y3=0.749182y(x3)=0.7491820.400000,y4=0.689680y(x4)=0.6896800.500000,y5=0.643470y(x5)=0.6434690.600000,y6=0.611189y(x6)=0.6111880.700000,y7=0.593415y(x7)=0.5934150.800000,y8=0.590672y(x8)=0.5906710.900000,y9=0.603431y(x9)=0.6034301.000000,y10=0.632122y(x10)=0.632121,写成向量的形式：,方程组和高阶方程的数值解法,各种方法都可以直接运用过来。以两个方程的方程组为例,Euler公式,Runge-Kutta公式向量形式,Runge-Kutta公式一般形式：,例题：求微分方程组,满足初始条件,在区间0,1上的数值解,取h=0.1.,容易计算该方程组的精确解为：y=sinxcosx+exz=sinx+cosx+ex1,x0=0.000000,y0=0.000000z0=1.000000 x1=0.100000,y1=0.205342,yy1=0.210000,z1=1.210175,zz1=1.200009x2=0.200000,y2=0.422806,yy2=0.440006,z2=1.441955,zz2=1.400139x3=0.300000,y3=0.654718,yy3=0.690043,z3=1.697616,zz3=1.600716x4=0.400000,y4=0.903649,yy4=0.960182,z4=1.979671,zz4=1.802304x5=0.500000,y5=1.172443,yy5=1.250564,z5=2.290899,zz5=2.005729x6=0.600000,y6=1.464238,yy6=1.561426,z6=2.634373,zz6=2.212097x7=0.700000,y7=1.782506,yy7=1.893128,z7=3.013488,zz7=2.422813x8=0.800000,y8=2.131082,yy8=2.246191,z8=3.431996,zz8=2.639604x9=0.900000,y9=2.514206,yy9=2.621320,z9=3.894045,zz9=2.864540 x10=1.000,y10=2.936564,yy10=3.019451,z10=4.404220,zz10=3.100055,经计算得数值解：,1、,2、确定方法，然后求解,（0.202760.0881157）（0.2130070.0934037）（0.2237630.0988499）（0.2350520.104437）（0.2469020.110146）,4阶Runge-Kutta法，h=1,高阶方程,则有：,令,用若干节点处的y及y值的线性组合来近似y(xn+1)。其通式可写为：,当10时，为隐式公式;1=0则为显式公式。,线性多步法,（1）求出开头几个点上的近似值，即计算“表头”；,线性多步法：,（2）利用逐步求后面点xk上的值yk。,基于数值积分的构造法,将在上积分，,只要近似地算出右边的积分，则可通过近似y(xn+1)。而选用不同近似式Ik，可得到不同的计算公式。,得到,积分系数,这是显格式，q+1阶r+1步格式。r=maxp,q,局部截断误差,例：建立p=1,q=2的显格式,p=1，,q=2，显格式，,积分区间为,积分节点为,所以,例：建立p=2,q=2的隐格式,p=2，,q=2，隐格式，,积分区间为,积分节点为,所以,它的截断误差较显格式小，通常也具有更好的稳定性。,阿当姆斯外推公式,以xn-2，xn-1，xn为节点作牛顿向后插值多项式P2(x)。,其中,数值解的稳定性,其中与f有关。若中不含yn+1，则方法是显式的，否则是隐式的，所以一般显式单步法表示为,为了叙述方便初值问题的单步法可以写成如下统一形式,例如，Euler方法中，有,对于不同的方法，计算值yn与准确解y(xn)的误差各不相同。所以有必要讨论各种方法的截断误差。我们称en=y(xn)-yn为某一方法在xn点的整体截断误差。显然，en不仅与xn这步的计算有关，它与以前各步的计算也有关，所以误差被称为整体的。分析和估计整体截断误差en是复杂的。为此，在前面我们给出了局部误差的概念。,定义对于任意固定的xn=x0+nh，若对于初值问题的显式单步法产生的近似解yn，均有yny(xn)(h0同时n)，则称该方法是收敛的。在定义中，xn是固定的点，当h0时有n，n不是固定的。显然，若方法是收敛的，则在固定点xn处的整体截断误差en=y(xn)-yn趋于零。下面给出方法收敛的条件。定理设初值问题的单步法是p阶的（p1），且函数满足对y的Lipschitz条件即存在常数L0，使,对一切成立，则方法收敛，且。,例：考察初值问题在区间0,0.5上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。,例：向后Euler公式的稳定性,误差方程：,考察隐式欧拉法,可见绝对稳定区域为：,注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。,3阶RungeKutta,显式14阶方法的绝对稳定区域为,插值公式的余项为,则积分公式的截断误差为,k=3时的外推公式为,余项为：,将差分表示成函数值的和的形式：,二阶阿当姆斯外推公式可改写为：,三阶阿当姆斯外推公式可改写为：,2.阿当姆斯内插公式,将被积函数用以xn-1，xn，xn+1为插值节点的内插多项式,得到：,k=1,k=2时,阿当姆斯外推法与内插联合起来,考虑常微分方程的边值问题：,其中p(x)，q(x)和f(x)均为a,b上给定的函数，,，为已知数。,假定p(x)、q(x)及f(x)均为a,b上充分光滑的函数，且q(x)0，这时，边值问题存在连续可微的解，且唯一。,解二阶常微分方程边值问题的差分法,用差分法解边值问题的主要步骤是：,（1）将区间a,b离散化；,（2）在这些节点上，将导数差商化，从而把微分方程化为差分方程；,(3)解差分方程实际上就是解线代数方程组。,将a,b区间用节点,分成N等分，其中x0=a与xN=b称为边界点，,而x1,x2,xN-1称为内点。,例9.7试用差分法解方程,解将0,1划分为四等分，即取，得五个节点,差分方程为,将它改写成,在每个内点列方程得,由追赶法公式解得：,y3=1.4855y2=1.2802y1=0.7753,对于一般的差分方程,仍然考虑最简单的模型，即只有初值产生误差，看看这个误差的传播。,差