常微分方程 5

常微分方程初值问题,1,一阶微分方程的初值问题,f(t,u)关于u满足Lipschitz条件,保证了方程组的初值问题有唯一解。,2,3,t0,t1,t2,tn,u(t),t,u,t1=t0+h,t2=t1+h,.ti+1=ti+h,.tn=tn-1+h,t1=t0+h,t2=t0+2*h,.ti=t0+i*h,.tn=t0+n*h,建立差分算法的两个基本的步骤：1.建立差分格式，对解的存在域剖分；采用不同的算法逼近误差截断误差（相容性）；数值解对真解的精度整体截断误差（收敛性）；数值解收敛于真解的速度；差分算法舍人误差（稳定性）.2.差分格式求解,4,差分方法的基本思想“就是以差商代替微商”,考虑如下两个Taylor公式：,（1）,（2）,从（1）得到：,5,从（2）得到：,从（1）-（2）得到：,从（1）+（2）得到：,6,欧拉法(Euler),0,t,u,T,1.在求解域上等距离分割：,7,2.在有：,3.应用时采用如下递推方式计算：,4.Euler法的几何意义,0,t,在递推的每一步，设定,过点,作的切线，该切线的,方程为：,即：,例对初值问题,用Euler法求解，用,即，,10,11,N=5;h=1/N;y=zeros(N,1);y(1)=1;fori=1:Ny(i+1)=y(i)+h*(2*(i-1)*h+y(i)endx1=(0:0.1:1)y1=3*exp(x1)-2*x1-2x=(0:h:1)plot(x,y,-*,x1,y1,-);,解析解为：,误差分析,构造算法后，这一算法在实际中是否可行呢？也就是说是否使计算机仿真而不失真，这还需要进一步分析。相容性收敛性稳定性,12,1.局部截断误差-相容性,为了分析分析数值方法的精确度，常常在,成立的假定下，,估计误差,这种误差称为“局部截断误差”，如图。,局部截断误差是以点,的精确解,为出发值，用数值方法推进到下一个点,而产生的误差。,13,14,定义：数值方法的局部截断误差为O(hp+1)，则称这种数值方法的阶数为p;如果p0,称常微分方程初值问题是相容的（Consistent），因为p是整数，所以p=1,15,收敛性研究的是误差累积产生的整体截断误差.收敛：对任意的t(t0,T，成立若此时，整体截断误差满足则称方法的收敛为p阶的,2.整体截断误差收敛性,整体截断误差是以点,的初始值,为出发值，用数值方法推进i+1步到点,，所得的近似值与精确值的偏差：,称为整体截断误差。,16,特例，,若不计初始误差，即,则,即,3.舍入误差稳定性,假设一个计算机仅表示4个数字（小数点后面），,那么,计算,19,3舍入误差稳定性,计算数值解的过程中,难免有舍入误差.稳定性就是讨论舍入误差是否会随着计算无限扩大地传递下去.数值方法稳定性指对初始误差的连续依赖性,我们的要求是：最初产生的小误差在以后的计算中虽然会传递下去，但不会无限制的扩大，这就是稳定性所描述的问题。下面引进稳定性的概念：,设由初值,得到精确解，,由初值,得到精确解，,若存在常数,和充分小的步长,使得,则称数值方法是稳定的,通常这里定义的稳定性指h0情况下的稳定性。,t,u,0,20,21,数值求解微分方程过程示意,微分方程,区域剖分,离散系统的性态研究,递推计算或解线性代数方程组,微分方程离散,初始和边界条件处理,解的存在性、唯一性,解的收敛性和收敛速度,解的稳定性,得到数值解,计算例题,其解析解为：,x=00.20000.40000.60000.80001.0000y=1.00001.20001.37331.53151.68111.8269,22,23,24,隐式欧拉公式假设使用向后差商替代方程中的导数项y(xn+1)，再离散化，即可导出下列格式用显式欧拉公式给出初始值，再由隐式公式进行迭代，因此得到如果迭代过程收敛，为隐式方程的解,25,隐式欧拉公式的局部截断误差,26,改进欧拉法,矩形法将方程的两端从xn到xn+1求积分得用左矩形法进行数值积分则因此有,27,梯形法用梯形公式计算上式右端的积分，可期望得到较高的精度，这时代入，从而有,28,2.3改进欧拉法,将上式中的y(xn)，y(xn+1)分别用yn，yn+1代替，作为离散化的结果导出下列计算公式将这个结果代入，并将其中的y(x1)用y1近似代替，则得一般地当求出yn以后，要求yn+1，则可归结为解方程,29,梯形公式的局部截断误差为,30,欧拉法与梯形法则相结合，计算公式为先用欧拉法由(xn,yn)得出y(xn+1)的初始近似值，然后进行迭代，反复改进这个近似值，直到(为所允许的误差)为止，并把取作y(xn+1)的近似值yn+1,31,在实际上，当h取值较小时，让梯形法的迭代公式只迭代一次就结束先用欧拉公式求得一个初步近似值，称之为预报值预报值的精度不高，用它替代梯形法右端的yn+1，再直接计算得出yn+1，并称之为校正值，这时得到预报校正公式,32,预报校正公式称为改进欧拉公式,这个公式还可写为,33,2.3改进欧拉法,34,2.3改进欧拉法,改进公式的截断误差,35,代入，可得y(xn+1)的二阶泰勒展开式为因此有,36,例：在区间0,1.5上，取h=0.1，求解,37,解：用欧拉法计算公式如下用迭代一次的改进欧拉法计算公式如下,38,2.3改进欧拉法,39,2.3改进欧拉法,例：用欧拉预报-校正公式求解初值问题要求步长，计算的近似值,40,2.3改进欧拉法,解：由题意知，改进欧拉法的具体形式为由，计算如下,41,2.3改进欧拉法,n=0时,42,2.3改进欧拉法,n=1时,三、改进的Euler法,将微分方程,在区间,上积分，得到,用梯形法计算积分的近似值，有,于是,这是一个隐式格式，一般需要用迭代法来求，而用显式的Euler法提供初值。,43,为了简化计算的过程，在此基础上进一步变为如下算法：,此式称为“改进的Euler法。,接下来讨论其几何意义,预估,校正,其局部截断误差为,这个问题将在下节讨论。,44,t,u,0,45,Euler法、改进的Euler法和解析解的比较,46,四、（龙格-库塔）Runge-Kutta方法,简单的Euler法是建立在Taylor级数的一项展开；,改进的Euler法是以两项Taylor级数为基础建立的，如：,如果我们截取Taylor级数的更多项会得到什么样的求解方法呢？,两个德国数学家（C.Runge&M.kutta）以这种思想为基础建立了求解微分,方程的龙格-库塔方法。它是常微分方程数值解法中使用最为广泛的方法之一。,47,用f(x,y)在几个不同点的数值加权平均代替，使截断误差的阶数尽可能高。也就是说，取不同点的斜率加权平均作为平均斜率，从而提高方法的阶数。,应用积分公式,无法计算,一般地，一个K阶的Runge-Kutta方法可用下面的公式表示：,其中，,是待定的加权系数，,是待定的系数。,Euler法就是,的R-K法。,其系数的确定如下：将,展开成,的幂级数，并与微分方程的精确解,在点,的Taylor展开式相比较，使两者的前,项相同，这样确定的R-K法，,其局部截断误差为,，根据所得关于待定系数的方程组，求出它们的值后,代入公式，就成为一个,阶R-K方法。,49,例题以二阶R-K法为例说明上述过程,把,代入,中，有,50,经比较得到,取为自由参数：,从而得到不同的但都是二阶的R-K方法，对应的有中点法、Heun（亨）法,以及改进的Euler法。,51,基于相同的过程，通过比较五次Taylor多项式，得到更加复杂的结果，给出了包含,13个未知数的11个方程。得到多组系数，其中常用的是以下四阶R-K法：,改进的Euler法、R-K法以及解析解的比较：,52,53,54,龙格库塔四阶法的特点龙格库塔四阶公式精度较高，可满足一般工程计算的要求每次计算yn+1时，只用到前一步的计算结果yn，因此在已知y0的条件下，可以自动的进行计算可以在计算过程中随时改变步长h缺点是每前进一步需要多次调用函数f(x,y)，工作量较大，并且误差不容易估计,55,7.3.5四阶龙格库塔法,56,3.5四阶龙格库塔法,例：写出四阶龙格-库塔法求解初值问题的计算公式，并取步长，计算的近似值解：由题意知，所以有,57,3.5四阶龙格库塔法,58,3.5四阶龙格库塔法,代入，有由于，取，取步长，有,两步欧拉公式为了提高精度，改用中心差商替代中的导数项，并离散化，有两步欧拉公式两步欧拉公式的局部截断误差,两步欧拉公式的局部截断误差是，是二阶方法,五、线性多步（LinearMultistepMethod）法,1.预备知识：插值多项式,插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。,从几何上理解：对一维而言，已知平面上n1个不同点，要寻找一条n次多项式曲线通过这些点。插值多项式一般常见的是拉格朗日插值多项式。,2.气象应用,不均匀站点上的气象要素数据,均匀网格点上的数据,插值,3.拉格朗日插值多项式,拉格朗日插值多项式逼近可能是求插值节点不均匀的插值多项式的最简单的方法。,实验观察结果或原始测量数据的分布通常是非均匀的。,例如，四个点可以确定一个三次多项式，其拉格朗日形式为：,61,4.Adams-Bashforth（阿达姆斯贝雪福斯）公式,首先，用以下四个点对,进行三次Langrage插值：,则,于是，有,容易算出，,例如，我们