常微分方程的数值解.ppt

常微分方程的数值解法,电子科技大学,常微分方程的数值解,引言 简单的数值方法 Runge-Kutta方法 一阶常微分方程组和高阶方程,在高等数学中我们见过以下常微分方程：,6.1 引言,(1)，(2)式称为初值问题，(3)式称为边值问题。,在实际应用中还经常需要求解常微分方程组：,本章主要研究问题(1)的数值解法，对(2)(4)只作简单介绍。,（其中L为Lipschitz常数）则初值问题（1）存在唯一的连续解。,考虑一阶常微分方程初值问题,其中，y = y(x) 是未知函数，y(x0) = y0 是初值条件，而f (x, y) 是给定的二元函数.,由常微分方程理论知，若f(x)在xa,b连续且 f 满足对 y 的Lipschitz条件：,常微分方程的数值解法有单步法和多步法之分： 单步法：在计算yn1 时只用到前一点yn 的值 ； 多步法：计算yn1 时不仅利用yn，还要利用yn-1, yn-2,. 一般k步法要用到 yn, yn-1, yn-2,., yn-k+1。,求问题（1）的数值解，就是要寻找解函数在一系列离散节点x1 x2 xn xn+1 上的近似值y1, y 2,yn 。,为了计算方便，可取xn=x0+nh，(n=0,1,2,)， h称为步长。,6.2 简单的数值方法,一、欧拉（Euler）方法,在x= x0 处，用差商代替导数：,由,得,同理，在x= xn 处，用差商代替导数：,由,得,若记,则上式可记为,此即为求解初值问题的Euler方法，又称显式Euler方法。,Euler方法的几何意义：,（Euler折线法）,例: 用Euler方法求解常微分方程初值问题,并将数值解和该问题的解析解比较。,解：Euler方法的具体格式：,xn y(xn) yn yn-y(xn) 0.0 0 0 0 0.2 0.1923 0.2000 0.0077 0.4 0.3448 0.3840 0.0392 0.6 0.4412 0.5170 0.0758 0.8 0.4878 0.5824 0.0946 1.0 0.5000 0.5924 0.0924 1.2 0.4918 0.5705 0.0787 1.4 0.4730 0.5354 0.0624,取h=0.2, xn=nh,(n=0,1,2,15), f(x,y)=y/x 2y2 计算中取f(0,0)=1. 计算结果如下：,xn y(xn) yn yn-y(xn) 1.6 0.4494 0.4972 0.0478 1.8 0.4245 0.4605 0.0359 2.0 0.4000 0.4268 0.0268 2.2 0.3767 0.3966 0.0199 2.4 0.3550 0.3698 0.0147 2.6 0.3351 0.3459 0.0108 2.8 0.3167 0.3246 0.0079 3.0 0.3000 0.3057 0.0057,由表中数据可以看到，微分方程初值问题的数值解和解析解的误差一般在小数点后第二位或第三位小数上，这说明Euler方法的精度是比较差的。, : 数值解 : 准确解,h=0.2;y(1)=0.2;x=h:h:3; for n=1:14 xn=x(n);yn=y(n); y(n+1)=yn+h*(yn/xn-2*yn*yn); end x0=0:h:3;y0=x0./(1+x0.2); plot(x0,y0,x,y,x,y,o),若直接对y=f(x,y)在xn, xn+1积分，,利用数值积分中的左矩形公式：,此即为Euler公式。,设y(xn)= yn，则得,若用右矩形公式：,得,上式称后退的Euler方法，又称隐式Euler方法。,可用迭代法求解：,初值：,迭代：,k=0,1,因,故当hL1时，迭代法收敛。,二、梯形方法,由,利用梯形求积公式：,得,上式称梯形方法，是一种隐式方法。,用迭代法求解：,初值：,迭代：,k=0,1,因,故当hL/21时，迭代法收敛。,由以上分析可以看出，隐式方法的计算比显式方法复杂，需要用迭代法求解非线性方程才能得出计算结果。,可采用将显式Euler格式与梯形格式结合使用的方法来避免求解非线性方程。,记,再用梯形格式计算：,预测,校正,上面两式统称预测校正法，又称改进的Euler方法。,三、单步法的局部截断误差和精度,单步法的一般形式为 ： （与f 有关）,显式单步法形式为：,整体截断误差：从x0开始，考虑每一步产生的误差，直到xn，则有误差,称为数值方法在节点xn处的整体截断误差。,但en不易分析和计算，故只考虑从xn到xn+1的局部情况。,定义：设y(x)是初值问题(1)的精确解，则称,为显式单步法在节点xn+1处的局部截断误差 。,若存在最大整数p使局部截断误差满足,则称显式单步法具有p阶精度或称p阶方法。,注：将Tn+1表达式各项在xn处作Taylor展开，可得具体表达式。,Euler方法的局部截断误差：,故Tn+1= O(h2)，p=1，,（设yn=y(xn)）,其中,称局部截断误差主项。,即Euler方法具1阶精度。,（设yn=y(xn)）,故Tn+1= O(h3)，p=2，,梯形方法的局部截断误差：,局部截断误差主项为：,梯形方法具2阶精度。,6.3 RungeKutta方法,一、Runge-Kutta方法的基本思想,由Taylor展式,Tn+1= O(hp+1)，若提高p，可提高精度。,但因,高阶导数计算复杂，故可从另外角度考虑。,分析Euler公式及改进的Euler公式：,局部截断误差：O(h2),局部截断误差：O(h3),可用f(x,y)在某些点处值的线性组合得yn+1，增加计算f(x,y)的次数可提高阶数。,设法计算f(x,y)在某些点上的函数值，然后对这些函数值作线性组合，构造近似计算公式，再把近似公式和解的泰勒展开式相比较，使前面的若干项吻合，从而获得达到一定精度的数值计算公式 。,Runge-Kutta方法的基本思想：,设,ci , i , ij 为待定常数。,上面第一个式子的右端在(xn,yn)作泰勒展开后,按h的幂次作升序排列 ：,再与初值问题的精确解y(x)在点x=xn处的泰勒展开式,相比较，使其有尽可能多的项重合。,例如，要求,就得到p个方程,从而定出参数ci ，i，ij ，再代入K1，K2， Kr的表达式，就可得到计算微分方程初值问题的数值计算公式 ：,若 Tn+1= O(hp+1)，则称其为p 阶r 级RK方法。,上式称为r 级RungeKutta方法的计算公式。,当r=1时，就是Euler方法。,要使Runge-Kutta公式具有更高的阶p，就要增加r 的值。下面我们只就r =2推导R-K方法。,二、二阶Runge-Kutta方法,其中 c1, c2, 2, 21 待定。,上式的局部截断误差为：,又,由,利用二元函数的Taylor展开，得,代入Tn+1的表达式，得,即,c1 = 1-a ， 2 = 21 =1/(2a),要使上式p=2阶，则需,方程组解不唯一，可令c2=a 0 ，则,满足上述条件的公式都为2阶R-K公式。,称中点公式，相当于数值积分的中矩形公式：,如取a= ，则c1= c2= ， 2=21 =1，即为改进Euler公式。,若取a= 1，则c1= 0，c2= 1， 2=21 = ，得,例:蛇形曲线的初值问题,令f(x,y)=y/x 2y2, 取f(0,0)=1, h=0.2, xn=nh , ( n = 1,2,15),2阶龙格-库塔公式计算格式： k1=yn/xn 2yn2, k2 = (yn+hk1)/(xn+h) 2(yn+hk1)2 yn+1=yn + 0.5h k1 + k2,x0=0;y0=0;h=.2; x=.2:h:3; k1=1; k2=(y0+h*k1)/x(1)-2*(y0+h*k1)2; y(1)=y0+.5*h*(k1+k2); for n=1:14 k1=y(n)/x(n)-2*y(n)2; k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)2; y(n+1)=y(n)+0.5*h*(k1+k2); end y1=x./(1+x.2); subplot(221) plot(x,y1,x,y1,b*) subplot(222) plot(x,y,x,y,o) subplot(223) plot(x,y,x,y,o,x,y1,x,y1,b*),三、三阶与四阶Runge-Kutta方法,当r=3时，R-K公式表示为,其中,为8个待定常数。,上式的局部截断误差为,类似二阶的推导过程，将K2, K3按二元函数展开，使Tn+1= O(h4)，得,方程有8个未知数，解不唯一。,满足该条件的公式统称为三阶R-K公式。,其中一个常用公式为：,当r=4时，利用相同的推导过程，经过较复杂的计算，可以得出四阶R-K公式的成立条件。,下列经典公式是其中常用的一个：,四、一阶常微分方程组和高阶微分方程的 数值解法简介,一阶常微分方程组的数值解法：,下列包含多个一阶常微分方程的初值问题：,称为一阶常微分方程组的初值问题。,引进向量记号:,则上述一阶常微分方程组的初值问题化为矩阵形式：,它在形式上跟单个微分方程的初值问题形式完全相同，只是函数变成了向量函数。故前面介绍的一切数值方法都适用，只要把函数换成向量函数即可。,k1=f(xn,yn), k2=f(xn+0.5h,yn+0.5hk1) k3=f(xn+0.5h,yn+0.5hk2), k4=f(xn+h,yn+hk3) yn+1= yn+ hk1+2k2+2k3+k4/6 (n = 0, 1, 2, ,N),四阶龙格-库塔公式,数值求解:,狐狸和野兔问题常微分方程组,在一个生物圈中,只有野兔和狐狸两种动物,记y1 为野兔数量, y2 为狐狸数量. 设有足够的食物供野兔生存,而狐狸只以野兔为食物. 模型如下,自变量x (0, 15), 初值条件为: y1(0)=20, y2(0)=20,t0=0;t1=15; %设置区域 y0=20 20; %给定初值 t,y=ode23(lot1,t0,t1,y0); %求解 plot(t,y) %绘图,function z=lot1(x,y) z(1,:)=y(1)-0.01*y(1).*y(2); z(2,:)=-y(2)+0.02*y(1).*y(2);,求解常微分方程: (1)定义函数; (2)调用ode23,一、高阶常微分方程的数值解法：,可化为一阶常微分方程组求解。,例如，二阶常微分方程初值问题,引进新的变量,令z=y，即可将上述二阶方程化为如下的一阶方程组的初值问题:,例 求下列高阶微分方程的数值解：,解：显然,假设,则,即二阶问题化为微分方程组的初值问题。,五、二阶常微分方程边值问题数值方法,考虑方程：,结合下述三种边界条件之一：,边界问题的解法：打靶法、有限差分法,(5.4)式中,它们分别称为第一、第二、第三边界问题。,打靶法,基本思路：将边值问题转化为初值问题考虑。或者说适当选择初始值使初值问题的解满足边值条件。然后用求解初值问题的任一种有效的数值方法求解。,以第一边界条件为例,考虑边值问题：,取y0=a , 考虑初值问题,待定，由数值解法求解(5.5)得到在 处的解 ，,使 ，这里 为给定允许误差界，就停止迭代改进, 输出作为数值解。,求解非线性方程,若 , 则得所求数值解。当,该方程可用二分法、正割法或Newton法等来求解。若求得,对第二类边界问题，也可转化为考虑初值问题(5.5)，取,对第三类边界问题，仍可转化为考虑初值问题(5.5)，取,以 为待定参数。,，以 为待定参数。,有限差分法,则离散化成差分方程,将区间a,b进行等分：,设在,处的数值解为 。,用中心差分近似微分，即,二阶中心差商,对应的边界条件也离散成,第