第三章常微分方程数值解课件.ppt

化工应用数学第三章,第三章,常微分方程的数值解法,主,要,内,容,3.1,、引言,3.2,、初值问题,3.3,、边值问题,?,1,化工应用数学第三章,3.1,、引,言,1.,什么是微分方程,?,在化学工程中关于扩散、反应、传质、传热,和流体流动等问题的数学模型。,微分方程：,是包含一个未知函数及其导数关,系的方程。,常微分方程：,其中只包含一个自变量的导数,的方程。,?,2,化工应用数学第三章,例题：,r,dc,p,P,?,dt,?,k,1,c,A,r,dc,s,A,?,k,S,?,1,?,P,dt,?,k,2,c,A,r,?,dc,T,T,A,?,k,2,?,S,dt,?,k,3,c,A,r,?,dc,A,A,?,k,A,?,dt,?,?,k,1,?,k,2,?,k,3,?,c,A,3,?,T,对于,r,?,dc,A,A,?,dt,?,?,k,1,?,k,2,?,k,3,?,c,A,:,求解,c,A,利用以前所学知识，进行积分可得：,c,(,k,A,?,e,?,1,?,k,2,?,k,3,),t,c,A,0,求解这样的方程一般我们需要给定初值。,?,3,化工应用数学第三章,例题：,一维均匀介质稳态导热问题。设其一端绝热，另一端恒,温为,T,1,，则此问题模型可用一常微分方程边值问题来描述。,d,dT,?,(,k,),?,0,?,?,dx,dx,?,dT,?,x,?,0,?,0,T,(,L,),?,T,1,?,?,dx,由以上两例可以看出，初值问题和边值问题区别在于：前者,在自变量一端给定附加条件，后者在自变量两端附加条件。,?,4,?,3.2,、初值问题,1.,尤拉法（,Euler,）,2,局部截断误差,3,改进尤拉法,4.,龙格,-,库塔法,5.,常微分方程组,6.,步长的选择,7.,收敛性和稳定性,8.,线性多步法,化工应用数学第三章,5,化工应用数学第三章,1.,尤拉,(Euler),法,一般初值问题可表示为,各,?,a,?,x,?,b,(1.1),种,?,y,?,?,f,(,x,y,),数,?,y,(,x,0,),?,y,值,0,(1.2),解,法,其中,f,(,x,y,),是已知函数，,(1.2),是定解条件也称为,初值,条件。,?,6,化工应用数学第三章,1.,尤拉,(Euler),法,考虑一阶初值问题,：,?,?,f,(,x,y,),y,?,?,?,y,(,a,),?,y,?,0,?,a,?,x,?,b,弄清常微方程初值,问题数值解法的一,些基本概念和构造,方法的思路,.,通过欧拉方法的讨论,最简单而直观实用方法,?,7,化工应用数学第三章,1.,尤拉,(Euler),法,把区间,a,b,分为,n,个小区间,步长为,h,i,?,(,x,i,?,1,-,x,i,),n,等份,节点,x,i,?,a,?,ih,i,一般取,h,i,?,h,(,?,(,b,?,a,),/,n,),即等距，要计算出解函数,y,(,x,),在一系列节点,a,?,x,0,?,x,1,?,.,?,x,n,?,b,处的近似值,u,0,u,1,?,u,n,?,8,化工应用数学第三章,1.,尤拉,(Euler),法,?,?,y,?,?,f,(,x,y,),a,?,x,?,b,?,y,(,a,),?,y,0,假定是,y(x),初值问题的解，那么，把在节点,x,i,附近展开,成,Taylor,级数，则,2,y,(,x,?,(,x,h,i,?,1,),?,y,(,x,i,),?,h,y,i,),?,2,y,?,?,(,?,),2,?,y,(,x,h,i,),?,hf,(,x,i,y,(,x,i,),?,2,y,?,?,(,?,),取前两项，并令,u,i,?,y,(,x,i,),u,i,?,1,?,u,i,?,hf,(,x,i,u,i,),?,9,化工应用数学第三章,1.,尤拉,(Euler),法,u,i,?,1,?,u,i,?,hf,(,x,i,u,i,),依上述公式逐次计算可得：,u,1,?,u,0,?,hf,(,x,0,u,0,),u,2,?,u,1,?,hf,(,x,1,u,1,),. . .,u,n,?,1,?,u,n,?,hf,(,x,n,u,n,),每步计算,u,n,?,1,只用到,u,n,?,10,化工应用数学第三章,几何意义,过点,(,x,0,y,0,),的曲线是解,y,(,x,),在,(,x,0,y,0,),作,y,(,x,),的切线,(,斜率,u,1,),与直线,x,?,x,1,交于,P,1,点,?,再作切线（斜率,u,2,）与直线,x,?,x,2,交于,P,2,点,?,?,?,P,2,P,1,y,?,y,(,x,),P,0,x,0,x,1,x,2,?,11,化工应用数学第三章,1.,尤拉,(Euler),法,例题：,用尤拉法解下列方程组，取,h,=0.1,。,?,?,?,y,?,?,y,?,2,x,y,(,0,?,x,?,1,),?,?,y,(,0,),?,1,.,解,取步长,h,=0.1,，欧拉公式的具体形式为,y,y,2,x,n,n,?,1,?,n,?,h,(,y,n,?,y,),n,已知,y,0,=1,由此式可得,y,?,h,(,y,2,x,0,1,?,y,0,0,?,y,),?,1,?,0,.,1,?,1,.,1,0,y,?,y,2,x,1,0,.,2,2,1,?,h,(,y,1,?,y,),?,1,.,1,?,0,.,1,(,1,.,1,?,.,1,),?,1,.,191818,1,1,?,.,12,化工应用数学第三章,1.,尤拉,(Euler),法,x,n,尤拉公式数值解,y,n,0.2,0.4,0.6,0.8,1.0,1.191818,1.358213,1.508966,1.649783,1.784770,准确解,y,(,x,n,),1.183216,1.341641,1.483240,1.612452,1.732051,误差,0.008602,0.016572,0.025726,0.037331,0.052719,与准确解,y,?,1,?,2,x,相比，可看出欧拉公式的计算结,果误差较大,.,?,13,化工应用数学第三章,2.,局部截断误差,u,i,?,y,(,x,i,),即第,i,步计算是精确的前提下，考虑的,定义,在假设,R,i,?,y,(,x,i,?,1,),?,u,i,?,1,称为局部截断误差。,截断误差,定义,若某算法的局部截断误差为,O,(,h,p,+1,),，则称该算法有,p,阶精度。,?,尤,拉法的局部截断误差：,R,i,?,y,(,x,i,?,1,),?,u,i,?,1,?,y,(,x,i,),?,hf,(,x,i,y,i,),?,h,2,2,y,?,?,(,x,i,),?,y,(,x,i,),?,hf,(,x,i,y,i,),2,?,h,2,2,y,?,?,(,x,i,),?,O,(,h,),因此，尤拉法具有,1,阶精度。,?,14,化工应用数学第三章,3.,改进的尤拉方法,显式尤拉公式,将,y,(,x,i,+1,) =,y,(,x,i,+,h,),在,处,Taylor,展开,隐式尤拉公式,将,y,(,x,i,+1,) =,y,(,x,i,+,h,),在,处,Taylor,展开,显式法计算较简单，隐式法计算复杂，可以编程。,?,15,化工应用数学第三章,3.,改进的尤拉方法,隐式尤拉公式,显式,与,隐式,两类方法各有特点，考虑到数值稳,定性等其他因素，人们有时需要选用,隐式,方法，但,使用,显式,算法远比,隐式,方便。,隐式方程通常用迭代法求解，而迭代过程的实,质是,逐步显式化,。,?,16,化工应用数学第三章,3.,改进的尤拉方法,尤拉法,隐式尤拉法,两式相加可得改进的尤拉公式，通常将其与显示尤拉公式联用,h,?,?,u,i,?,1,?,u,i,?,f,(,x,i,y,i,),?,f,(,x,i,?,1,y,i,?,1,),?,?,2,?,u,i,?,1,?,u,i,?,hf,(,x,i,y,i,),?,?,17,化工应用数学第三章,3.,改进的尤拉方法,?,改进尤拉法,/* modified Eulers method */,Step 1,:,先用,显式尤,拉公式作,预测,，算出,u,i,?,1,?,u,i,?,hf,(,x,i,y,i,),y,Step 2,:,再将,1,代入,隐式,梯形公式的右边作,校正,，得到,i,?,u,i,?,1,h,?,u,i,?,f,(,x,i,y,i,),?,f,(,x,i,?,1,y,i,?,1,),2,h,u,i,?,1,?,u,i,?,?,f,(,x,i,y,i,),?,f,?,x,i,?,1,u,i,?,h,f,(,x,i,y,i,),?,?,2,(,i,?,0,.,n,?,1,),改进的尤拉公式是二阶方法，但计算量增加了一倍。,?,18,化工应用数学第三章,例：,?,?,y,?,x,?,y,?,1,x,?,0,0.5,?,y,(0),?,1,h,?,0.1,解：由改进的尤拉公式可得,u,?,u,h,i,?,1,i,?,2,f,(,x,i,y,i,),?,f,(,x,i,?,1,y,i,?,1,),y,i,?,1,?,y,?,0.1,i,2,x,i,?,y,i,?,1,?,x,i,?,1,?,y,i,?,0.1(,x,i,?,y,i,将数值代入，进行逐步计算。,一般对于改进的尤拉法可采用计算机编程进行计算。,?,?,1),?,1,19,化工应用数学第三章,4.,龙格,-,库塔法,建立高精度的单步递推格式。,基本思想,:,是从,(,x,i,y,i,),点出发，以,某一斜,率,沿直线,达到,(,x,i,+1,y,i,+1,),点。欧拉法及其各种变形所能达到的最高精度为,2,阶。,?,考察改进的尤拉法，可以将其改写为：,y,i,?,1,K,1,K,2,?,?,?,?,1,?,1,?,y,i,?,h,?,K,1,?,K,2,?,2,?,2,?,f,(,x,i,y,i,),f,(,x,i,?,h,y,i,?,hK,1,),20,化工应用数学第三章,K,2,?,?,m,K,m,y,i,?,1,?,h,?,1,K,1,?,?,2,?,.,?,y,i,K,1,?,f,(,x,i,y,i,),?,2,h,y,i,?,?,21,hK,1,),K,2,?,f,(,x,i,?,?,3,h,y,i,?,?,31,hK,1,?,?,32,hK,2,),K,3,?,f,(,x,i,?,.,.,?,m,h,y,?,?,m,1,2,?,?,m m,?,1,?,K,m,?,f,(,x,i,?,?,?,m,?,.,),hK,m,1,1,hK,hK,2,其中,?,i,(,i,= 1, ,m,),，,?,i,(,i,= 2, ,m,),和,?,ij,(,i,= 2, ,m,;,j,= 1, ,i,?,1,),均为待定,系数，确定这些系数的,步骤与前面相似。,?,最常用为四级,4,阶,经典龙格,-,库塔法,/* Classical Runge-Kutta,Method */,：,y,i,?,1,K,1,K,2,K,3,K,4,?,?,?,?,?,y,i,?,h,(,K,1,?,2,K,2,?,2,K,3,?,K,4,),6,f,(,x,i,y,i,),h,f,(,x,i,?,h,y,?,K,1,),i,2,2,h,f,(,x,i,?,h,y,?,K,2,),i,2,2,f,(,x,i,?,h,y,i,?,hK,3,),?,23,化工应用数学第三章,x,?,0,0.5,?,y,?,x,?,y,?,1,例题,?,y,(0),?,1,h,?,0.1,解：,?,?,y,i,?,1,?,y,i,?,h,6,(,K,1,?,2,K,2,?,2,K,3,?,K,4,),K,1,?,f,(,x,i,y,i,),?,x,i,?,y,i,?,1,K,h,h,h,h,2,?,f,(,x,i,?,2,y,i,?,2,K,1,),?,(,x,i,?,2,),?,(,y,i,?,2,K,1,),?,1,K,h,h,h,h,3,?,f,(,x,i,?,2,y,i,?,2,K,2,),?,(,x,i,?,2,),?,(,y,i,?,2,K,2,),?,1,K,4,?,f,(,x,i,?,h,y,i,?,hK,3,),?,(,x,i,?,h,),?,(,y,i,?,hK,3,),?,1,K,1,?,0,K,2,?,0.05,K,3,?,0.0475,K,4,?,0.09525,将,K,值代入式，即可计算出方程组的解。,24,化工应用数学第三章,注：,?,龙格,-,库塔法的主要运算在于计算,K,i,的值，即计算,f,的,值。,Butcher,于,1965,年给出了计算量与可达到的最高精,度阶数的关系：,每步须算,K,i,的个数,2,3,4,5,6,7,n,?,8,可达到的最高精度,O,(,h,2,),O,(,h,3,),O,(,h,4,),O,(,h,4,),O,(,h,5,),O,(,h,6,),O,(,h,n,?,2,),?,25,化工应用数学第三章,5.,常微分方程组,?,?,y,(,x,),?,f,(,x,y,z,),，,y,(,x,0,),?,y,0,?,z,(,x,),?,g,(,x,y,z,),，,z,(,x,0,),?,z,0,?,?,?,y,y,h,i,?,1,?,i,?,(,k,1,?,2,k,2,?,2,k,3,?,k,4,),?,6,?,?,?,z,i,?,1,?,z,i,?,h,6,(,L,1,?,2,L,2,?,2,L,3,?,L,4,),?,?,?,k,1,?,f,(,x,i,y,i,z,i,),?,?,k,h,h,h,2,?,f,(,x,i,?,L,z,?,?,2,y,i,?,2,k,1,z,i,?,2,L,1,),1,?,g,(,x,i,y,i,i,),?,?,?,L,2,?,g,(,x,i,?,h,2,y,i,?,h,2,k,1,z,i,?,h,2,L,1,),?,?,?,k,?,f,(,x,h,h,k,h,3,i,?,y,i,?,2,z,i,?,L,2,),?,2,2,2,?,?,?,k,4,?,f,(,x,i,?,h,y,i,?,hk,3,z,i,?,hL,3,),?,?,L,3,?,g,(,x,i,?,h,2,y,i,?,h,2,k,2,z,i,?,h,2,L,2,),?,L,4,?,g,(,x,i,?,h,y,i,?,hk,3,z,i,?,hL,3,),?,26,化工应用数学第三章,5,.,常微分方程组,高阶常微分方程,?,?,y,(,x,),?,f,(,x,y,y,),?,y,(,x,?,0,),?,y,0,?,y,(,x,0,),?,y,1,思路：转化降阶，通常转化为常微分方程组（,令,z=y,?,?,?,z,(,x,),?,f,(,x,y,z,),z,(,x,0,)=,y,1,?,y,(,x,),?,z,y,(,x,0,),?,y,0,?,1,阶）来求解。,27,化工应用数学第三章,?,y,?,2,y,?,2,y,?,e,2,x,sin,x,例题,?,?,y,(0),?,?,0.4,?,?,y,(0),?,?,0.6,令,z,=,y,则方程转化为,?,z,?,2,z,?,2,y,?,e,2,x,?,sin,x,z,(0),?,?,0.6,?,z,=,y,y,(0),?,?,0.4,?,?,y,?,z,?,?,?,?,y,h,i,?,1,?,y,i,?,(,k,1,?,2,k,2,?,2,k,3,?,k,4,),?,?,z,?,2,z,?,2,y,?,e,2,x,sin,x,?,6,?,?,?,z,i,?,1,?,z,i,?,h,6,(,L,1,?,2,L,2,?,2,L,3,?,L,4,),?,28,化工应用数学第三章,?,?,?,k,1,?,z,i,?,L,2,x,1,?,e,i,sin,x,i,?,2,z,i,?,2,y,i,?,?,?,k,?,h,2,?,z,i,?,2,L,1,?,2(,x,h,i,?,?,?,L,2,?,e,2,),sin(,x,h,1,h,i,?,2,),?,2(,z,i,?,2,hL,1,),?,2(,y,i,?,2,k,1,),?,?,?,k,?,z,h,3,i,?,L,2,?,2,?,?,?,L,3,?,e,2(,x,h,i,?,2,),sin(,x,h,1,h,i,?,2,),?,2(,z,i,?,2,hL,2,),?,2(,y,i,?,2,k,2,),?,?,h,?,k,4,?,z,i,?,2,L,3,?,?,L,4,?,e,2(,x,i,?,h,),sin(,x,i,?,h,),?,2(,z,i,?,hL,3,),?,2(,y,i,?,hk,3,),29,化工应用数学第三章,6,.,步长的选择,从局部截断误差上来看，与,O,(,h,p,?,1,),有关，故,h,越小，对降低,局部误差有利，但是步长越小，势必引起计算量的加大。,利用,p,阶方法可构造出：,步长为,h,:,y,(,x,h,i,?,1,),?,y,i,?,1,?,ch,p,?,1,?,O,(,h,p,?,2,),（,1,）,步长为,h/2,:,h,y,(,x,i,?,1,),?,y,2,2,c,(,h,p,?,1,p,?,2,i,?,1,?,2,),?,O,(,h,),（,2,）,上述两式整理可得：,h,p,h,y,2,y,2,i,?,1,?,y,i,?,1,i,?,1,?,?,2,p,?,1,（,3,）,30,化工应用数学第三章,6,.,步长的选择,h,h,h,将,(3),式变形为：,y,(,x,1,i,?,1),?,y,2,y,2,i,?,1,?,y,i,?,i,?,1,?,2,p,?,1,h,令,?,=,y,2,i,?,1,?,y,h,i,?,1,(1),对于给定的精度,?,，如果,?,?,?,，可以将步长反复减半进行计算，,直到,?,?,?,为止，这时以最终得到的,y,i,?,1,作为结果；,(2),如果,?,?,?,，可以反复将步长加倍，直到,?,?,这时将步长再减,半一次，即为所求结果。,?,31,化工应用数学第三章,7.,收敛性和稳定性,收敛性与稳定性从两个不同的角度描述了微分方程数值解法,的实用价值，只有既收敛又稳定的方法，才可以提供比较可靠,的计算结果。,当步长,h,趋于零时，方程的数值解,y,i,?,1,趋近于,收敛性：,精确解,y,(,x,i,?,1,),时，则称所采用的数值方法是收敛的。,当步长,h,确定之后，随着步数的增加，计算中,稳定性：,累积的误差不会超出所允许的范围。,?,32,化工应用数学第三章,3.3,、边值问题,?,?,?,1.,打靶法,2,差分法,34,化工应用数学第三章,边值问题的数值解,/* Boundary-Value Problems */,2,阶常微分方程边值问题,x,?,(,a,b,),?,y,?,?,?,f,(,x,y,y,?,),?,?,y,(,a,),?,?,y,(,b,),?,?,?,打靶法,/* shooting method */,先猜测一个初始斜率,y,?,(,a,) =,m,，通过解初,值问题,y,?,?,?,f,(,x,y,y,?,),?,?,?,y,(,a,),?,a,?,y,?,(,a,),?,m,?,y,?,(,s,0,),斜率,=,s,0,y,(,x,),?,?,(,s,1,),y,(,b,) =,?,(,s,),找出,s,*,使得,?,(,s*,) =,?,，即把问,题转化为求方程,?,(,s,),?,?,= 0,的根。,?,斜率,=,s,1,0,a,b,x,35,1.,打靶法,例题：,?,1,?,?,y,?,8,?,y,?,4,?,y,(0),?,0,y,(10),?,0,用打靶法解。,?,化工应用数学第三章,36,化工应用数学第三章,1.,打靶法,?,?,?,?,y,?,8,?,1,y,?,4,?,y,(0),?,0,y,(0,）,=,m,3.744273,首先取初值,m,0,?,10,，利用四阶龙格,-,库塔法求出,y,0,(10),?,A,.,m,1,?,11,求出,y,1,(10),?,B,1.826444,然后利用,m,0,m,1,用弦截法求出下一个初值,m,2.,则,m,?,m,B,2,1,?,（,11,?,10,）求出,y,2,(10),?,1.08853,?,10,?,5,B,?,A,?,37,化工应用数学第三章,2,差分法,?,有限差分法,基本思想：,运用数值微分将导数用离散点上函数值表示，从而将边,值问题的微分方程和边界条件转化为只含有限个未知数的差分方,程组，并将此差分方程组的解作为该边值问题的数值解。,1.,二阶常微分方程的第一边值问题,?,?,?,y,?,(,x,),?,q,(,x,),y,(,x,),?,f,(,x,),a,?,x,?,b,?,y,(,a,),?,?,y,(,b,),?,?,其中,q(x)(,?,0),，,f(x),在,a,，,b),上连续，,?,，,?,为常数。,?,38,化工应用数学第三章,2,差分法,?,设等距节点：,x,i,= a + ih,，,i = 0,，,1,，,2,，,，,n,，,?,对其中内节点应用三点微分公式：,y,(,x,i,?,1,),?,2,y,(,x,i,),?,y,(,x,i,?,y,?,?,(,x,?,1,),2,i,),h,2,?,O,(,h,).,?,?,y,(,x,i,?,1,),?,2,y,(,x,i,),?,y,(,x,i,?,1,),2,h,2,?,q,(,x,i,),y,(,x,i,),?,f,(,x,i,),?,O,(,h,),?,(,当,h,充分小时，略去,O(h,2,),，并以,y,i-1,，,y,i,，,y,i+1,，代,y(x,i-1,),，,y(x,i,),，,y(x,i+1,),，得计算,y,i,的差分方程组,),?,y,i,?,1,?,2,y,i,?,y,i,?,1,?,h,2,?,q,i,y,i,?,f,i,39,化工应用数学第三章,2,差分法,加上边界条件即得边值问题的差分方程组,?,?,?,y,2,2,i,?,1,?,(,2,?,h,q,i,),y,i,?,y,i,?,1,?,h,f,i,i,?,1,2,?,?,?,n,?,1,?,y,0,?,?,y,n,?,?,其中,q,i,= q(x,i,),，,f,i,= f (x,i,),。,矩阵形式：,?,?,1,?,?,?,?,?,1,2,?,h,2,q,?,1,?,?,y,0,?,?,?,1,?,?,y,h,2,f,?,1,?,?,?,?,1,2,?,h,q,?,2,?,1,?,?,y,?,1,2,?,2,?,?,?,h,2,f,2,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,.,?,?,1,2,?,h,2,q,n,?,1,?,1,?,?,y,n,?,1,?,?,?,?,?,1,?,?,?,?,?,?,?,?,y,?,h,2,f,n,?,1,?,n,?,?,?,?,?,?,?,?,40,化工应用数学第三章,例,取,h = 0.1,，求边值问题：,?,?,y,?,?,(,x,),?,y,(,x,),?,x,0,?,x,?,1,的数值解。,?,y,(,0,),?,0,y,(,1,),?,1,解：本例的节点,x,i,= 0.1i,，,i = 0,，,1,，,2,，,，,10,，,q(x) = 1,，,r(x) = x,，,差分方程为,?,?,y,?,(2,?,h,2,),y,?,y,?,h,2,?,i,?,1,i,i,?,1,x,i,?,?,y,0,?,?,y,n,?,?,?,y,2,2,?,i,?,1,?,(,2,?,0,.,1,),y,i,?,y,i,?,1,?,0,.,1,x,i,i,?,1,2,?,?,?,n,?,1,?,y,0,?,0,y,n,?,1,?,42,化工应用数学第三章,差分方程为,?,?,y,i,?,1,?,2,.,01,y,i,?,y,i,?,1,?,0,.,01,i,i,?,1,2,?,?,?,n,?,1,?,y,0,?,0,y,n,?,1,?,?,?,2,.,01,1,?,?,y,?,1,?,2,.,01,1,?,?,1,?,?,?,?,0,.,001,?,y,?,?,?,2,?,.,.,.,?,?,?,?,?,?,?,0,.,002,?,.,?,?,?,1,?,2,.,01,1,?,?,?,?,y,?,?,?,n,?,2,?,?,?,0,.,008,?,?,1,?,2,.,01,?,?,?,?,y,n,?,1,?,?,?,?,?,0,.,991,?,?,解此方程组。,?,43,