数值分析课件ch09常微分方程数值解

1,第九章常微分方程数值解,2,常微分方程初始问题,许多科学技术问题，例如天文学中的星体运动，空间技术中的物体飞行，自动控制中的系统分析，力学中的振动，工程问题中的电路分析等，都可归结为常微分方程的初值问题。常微分方程中往往只有少数较简单和典型的常微分方程（线性常系数常微分方程等）可求出其解析解。本章考虑的是一阶微分方程（组）的数值解。高阶的常微分方程可以转化成一阶微分方程组。,3,举例,例：单摆运动问题,O,A,m,可以转化成一阶微分方程组,4,常微分方程解的存在唯一性,定理考虑一阶常微分方程的初值问题,对任意定义在上的都成立，则上述方程存在唯一解。,只要在上连续,且关于y满足Lipschitz条件，,即存在与无关的常数L使,5,常微分方程的数值解,数值解的基本思想,满足微分方程,常微分方程等价于下列积分方程,常微分方程初值问题的数值解法，是要寻求解函数y(x)在一系列节点处的近似值,关键在于导数如何近似。,关键在于积分如何近似。,6,初始问题的欧拉方法,初值问题的数值解法一般按节点从左至右的顺序依次求出的近似，所以称为步进法,单步法,从初值开始,依次求出,后一步的值只依靠前一步的来计算,典型的单步法是Euler(欧拉)方法,其计算格式是:,多步法在计算时需要前面若干步值来计算。,7,欧拉方法的推导,以Euler法为例说明初值问题数值方法的三种基本途径,Taylor展开法,忽略高阶项,取近似值可得到Euler公式,数值微分法,用差商代替导数,8,欧拉方法的推导,则,数值积分法,将区间积分,右端采用左矩形公式数值积分，得,9,几个简单的数值方法,显式欧拉法,隐式欧拉法,方法上是向后差商近似导数，或右矩形数值积分得到,由于未知数yn+1同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/欧拉公式.,一般先用显式计算一个初值，再迭代求解。,10,几个简单的数值方法,梯形公式,显、隐式两种算法的平均或梯形公式数值积分得到,改进的欧拉公式,11,单步法的误差分析,整体截断误差,数值解和精确解之差,整体截断误差除与步计算有关外,还与的计算有关，所以一般很难分析清楚。,局部截断误差,一般单步法有增量形式,在假设yn=y(xn)，即第n步计算是精确的前提下，考虑截断误差Tn+1=y(xn+1)yn+1称为局部截断误差,12,算法的精度,若某算法的局部截断误差为O(hp+1),则称该算法有p阶精度。,显式欧拉，隐式欧拉为1阶精度梯形公式和改进欧拉公式为2阶精度,显式欧拉局部截断误差,梯形公式局部截断误差,13,例,例：考察初值问题在区间0,0.5上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。,1.00002.00004.00008.00001.60001013.2000101,1.00002.50001016.25001021.56251023.90631039.7656104,1.00002.50006.25001.56261013.90631019.7656101,1.00004.97871022.47881031.23411046.14421063.0590107,从表中可以看出，这几种方法计算误差比较大，说明这些方法精度不够，需寻找更高精度的方法,14,高精度方法,基于Taylor展开式可以建立任意精度的单步法,而,所以，可以构造格式,这种格式使用到了各阶偏导数，使用不便。,15,Runge-Kutta方法,Runge-Kutta方法的基本思想,只要对平均斜率提供一种算法，便相应地得到一种微分方程的数值计算公式。,16,Runge-Kutta方法,考察改进的欧拉法,可以将其改写为：,改进欧拉公式却是利用了xn与xn-1两个点的斜率值k1=f(xn,yn)与k2=f(xn+1,yn+hk1)取算术平均作为平均斜率的近似值。,17,Runge-Kutta方法,启示设法在xk,xk+1上多预报几个点的斜率值，然后将它们加权平均作为平均斜率的近似值，则有可能构造出更高精度的计算公式,推广,其中由f(x,y)在一些点的值加权平均所得,这种方法称为Runge-Kutta方法,简记为R-K公式.,若由r个f(x,y)的值加权平均所得的公式称为r级R-K公式,18,二级Runge-Kutta方法,二级Runge-Kutta方法,希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得,将K2在(xi,yi)点作Taylor展开,19,二级Runge-Kutta方法,将K2代入第1式，得到,将yi+1与y(xi+1)在xi点的泰勒展开作比较,20,二级Runge-Kutta方法,要求，则必须有：,这里有3个未知数，2个方程,解不唯一,就是改进的欧拉法。,就是中点公式。,21,三阶与四阶显式R-K方法,三级Runge-Kutta方法,其中及均为待定参数.,局部截断误差为,只要将按二元函数泰勒展开，使，,22,三阶与四阶显式R-K方法,可得待定参数满足方程,这是8个未知数6个方程的方程组，解也不是惟一的.,23,三阶与四阶显式R-K方法,一个特殊的方法,同样的方法，可以得到经典4阶R-K方法,24,注,更一般R-K方法,龙格-库塔法的主要运算在于计算的值，即计算f的值。Butcher于1965年给出了计算量与可达到的最高精度阶数的关系：,25,注,由于龙格-库塔法的导出基于泰勒展开，故精,度主要受解函数的光滑性影响。对于光滑性不,太好的解，最好采用低阶算法而将步长h取小。,26,例,解,这里，经典的四阶龙格-库塔公式为,27,28,F=y-2*x/y;a=0;b=1;h=0.1;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=1;fori=1:nx=X(i);y=Y(i);K1=h*eval(F);,x=x+h/2;y=y+K1/2;K2=h*eval(F);x=x;y=Y(i)+K2/2;K3=h*eval(F);x=X(i)+h;y=Y(i)+K3;K4=h*eval(F);Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;end,29,30,求微分方程（组）的数值解命令,t，x=solver（f,ts,x0,options）,31,单步法的收敛性,定义若某算法对于任意固定的x=x0+nh，当h0(即n)时有yny(xn)，则称该算法是收敛的。,例,就初值问题考察欧拉显式格式的收敛性。,解：该问题的精确解为,欧拉公式为,对任意固定的x=xn=nh，有,32,定理:对初值问题的单步法,若局部截断误差为,且函数对y满足Lipschitz条件,即存在L0,使得,对一切成立,则该方法收敛,且有,单步法的收敛性,收敛性的条件,欧拉，改进欧拉以及Runge-Kutta方法都是收敛的整体截断误差比局部截断误差低一阶,33,对于欧拉方法，由于其增量函数就是，故当关于满足利普希茨条件时它是收敛的.,改进的欧拉方法，其增量函数由,欧拉方法的收敛性,改进欧拉方法的收敛性,假设关于满足利普希茨条件，记利普希茨常数为,34,单步法的稳定性,考虑收敛性时，当h0时有yny(xn)，此时yn为理论数值解，也就没有考虑到计算过程中的舍入误差当考虑计算过程中舍入误差对最后结果是否有影响，所以引入数值算法绝对稳定的概念,定义,由于f(x,y)比较复杂，给稳定性研究带来困难，所以考虑试验方程,35,单步法的稳定性,试验方程,该试验方程为线性方程，其精确解是稳定的。若对该方程数值算法都稳定，则对其它方程也靠不住。对一般非线性方程，局部可以化成线性的来近似,单步法的稳定性单步法用于试验方程，若得到的解为满足则是绝对稳定的。,36,单步法的稳定性,记，的范围叫绝对稳定区域,绝对稳定区域与实轴之交叫绝对稳定区间,方法A比方法B的绝对稳定区域大称A比B稳定,显式欧拉法的稳定性,在复平面上是以(-1,0)为圆心的单位圆,绝对稳定区间为【-2,0】,37,隐式欧拉法的稳定性,所以绝对稳定区域为左半平面绝对稳定区间隐式欧拉比显式欧拉稳定性好,38,梯形公式的稳定性,所以绝对稳定区域为左半平面绝对稳定区间,39,经典4阶龙格-库塔方法的稳定性,40,研究一般方程的稳定性时，相当于模型方程中的,41,线性多步法,一个特例,对于初值问题,若用中心差商离散导数,称为二步法,42,线性多步法,推广,当10时，为隐式公式;1=0则为显式公式。,43,线性多步方法的构造,基于Taylor展开的构造方法:,系数及可根据方法的局部截断误差及阶确定，其定义如下.,44,线性多步方法的构造,45,线性多步方法的构造,要使线性多步法达到p阶精度，则满足,上式为线性方程组，可能解不唯一。,46,一些特殊的线性多步法,Adams4步4阶显式公式,Adams3步4阶隐式公式,47,一些特殊的线性多步法,Milne4步4阶显