最新常微分方程的数值解PPT课件.ppt

常微分方程的数值解常微分方程的数值解常微分方程的数值解常微分方程的数值解n引言引言n简单的数值方法简单的数值方法nRunge-Kutta方法方法n一阶常微分方程组和高阶方程一阶常微分方程组和高阶方程例例:用用Euler方法求解常微分方程初值问题方法求解常微分方程初值问题并将数值解和该问题的解析解比较。并将数值解和该问题的解析解比较。解：解：Euler方法的具体格式：方法的具体格式：xn y(xn)yn yn-y(xn)0.00000.20.19230.20000.00770.40.34480.38400.03920.60.44120.51700.07580.80.48780.58240.09461.00.50000.59240.09241.20.49180.57050.07871.40.47300.53540.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.60.44940.49720.04781.80.42450.46050.03592.00.40000.42680.02682.20.37670.39660.01992.40.35500.36980.01472.60.33510.34590.01082.80.31670.32460.00793.00.30000.30570.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);endx0=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 k=0,1,=0,1,因因故当故当hL1时，迭代法收敛。时，迭代法收敛。二、梯形方法二、梯形方法由由利用梯形求积公式：利用梯形求积公式：得得上式称上式称梯形方法梯形方法，是一种隐式方法。，是一种隐式方法。用迭代法求解：用迭代法求解：初值：初值：迭代：迭代：k k=0,1,=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+k2x0=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);endy1=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)=20t0=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 a,考考虑虑初初值问题值问题 待定，由数值解法求解待定，由数值解法求解(5.5)得到在得到在 处的处的解解 ，使使 ，这里，这里 为给定允许误差界，就为给定允许误差界，就停止迭代改进停止迭代改进,输出作为数值解。输出作为数值解。求解非线性方程求解非线性方程若若 ,则得所求数值解。当则得所求数值解。当该方程可用二分法、正割法或该方程可用二分法、正割法或Newton法等来求解。法等来求解。若求得若求得对第二类边界问题，也可转化为考虑初值问题对第二类边界问题，也可转化为考虑初值问题(5.5)，取，取对第三类边界问题，仍可转化为考虑初值问题对第三类边界问题，仍可转化为考虑初值问题(5.5)，取，取 以以 为待定参数。为待定参数。，以，以 为待定参数。为待定参数。有限差分法有限差分法则离散化成差分方程则离散化成差分方程将区间将区间a,b进行等分：进行等分：设在设在处的数值解为处的数值解为 。用中心差分近似微分，即用中心差分近似微分，即 二阶中心差商二阶中心差商 对应的边界条件也离散成对应的边界条件也离散成 第一边界问题：第一边界问题：第二边界问题：第二边界问题：第三边界问题：第三边界问题：其中其中 为已知函数，则由常微分方程的理论知，通过为已知函数，则由常微分方程的理论知，通过则近似差分方程成离散差分方程为则近似差分方程成离散差分方程为其中其中变量替换总可以消去方程中的变量替换总可以消去方程中的 项，不妨设变换后的方程为项，不妨设变换后的方程为若若 是是 的线性函数时，的线性函数时，f 可写成可写成将以上方程合并同类项整理得方程组将以上方程合并同类项整理得方程组其中只要其中只要 ，则方程组的系数矩阵为，则方程组的系数矩阵为弱对角占优的三对角阵弱对角占优的三对角阵，其中其中 。而且还有而且还有误差估误差估 若若 不是不是 的线性函数时，所得方程组是非线性的线性函数时，所得方程组是非线性方程组为三对角线方程组，可以用方程组为三对角线方程组，可以用追赶法追赶法求解。求解。计计：组，可以用解非线性方程组的方法求解。例如，可用组，可以用解非线性方程组的方法求解。例如，可用简单迭代法简单迭代法、Newton迭代法迭代法求解。求解。