1、第12次 常微分方程初值问题数值解法计算方法计算方法(Numerical Analysis)1专业课件内容1.常微分方程初值问题解的存在性定理2.Euler公式3.梯形公式4.两步Euler公式5.欧拉法的局部截断误差6.改进型Euler公式7.龙格-库塔法8.算法实现2专业课件常微分方程初值问题解的存在性定理3专业课件第9章 常微分方程初值问题数值解法 包含自变量、未知函数及未知函数的导数的方程称为微分方程。 微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。(n)y,y,y 都是一次的,则称其为线性的,否则称为非线性的。 如果未知函数y及其各阶导数9.1 引言 自变量个数只有一个
2、的微分方程称为常微分方程。4专业课件 如下是一些典型方程求解析解的基本方法 可分离变量法、 常系数齐次线性方程的解法、 常系数非齐次线性方程的解法等。22yxy的解就不能用初等函数及其积分来表达。 但能求解的常微分方程仍然是很少的，大多数的常微分方程是不可能给出解析解。例如，一阶微分方程5专业课件 从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。00y)y(xy)f(x,y (9.1 ) 在区间a x b上的数值解法。 本章主要讨论一阶常微分方程初值问题 6专业课件定理1：如果函数f(x, y)在带形区域Ry,y,yyL)yf(x,)yf(x,212121则方程( 9.1 ) 在
3、 a, b 上存在唯一的连续可微分的解的解 y=y(x) 。 yb,-xaR内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与x, y无关)使 7专业课件推论：如果函数f(x,y)对y的偏导数 在带形区域L|yy)f(x,|对R内的所有x，y 都成立。yy)f(x,即存在常数L(它与x,y无关)使 则方程( 9.1 ) 在 a , b 上存在唯一的连续可微解y=y(x) 。yb,-xaR内有界。Home8专业课件Euler公式本章假设微分方程初值问题(9.1)有解9专业课件常微分方程初值问题(9.1)的数值解法的基本思想：算出精确解y(x)在区间 a,b 上的一系列离散
4、节点)y(x,),y(x),y(xn10n10y,y,y的近似值处的函数值 y=y(x)a=x0 xn=bx1x2x3(未知) bxxxxan1n1010专业课件i1ixxh 相邻两个节点的间距 称为步长，步长可以相等，也可以不等。数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。 a=x0 xn=bx1x2x3xn-1x4n,1,2,i ih,xx0i 本章总是假定h为定数，称为定步长，这时节点可表示为11专业课件 常微分方程数值解法的基本出发点：离散化。采用“步进式”，即求解过程顺着节点排列的次序逐步向前推进。y00y)y(xy)f(x,y中的导数 进行离散化处理。 以便对初
5、值问题0 02 2i i1 1i ii i, ,y y, , ,y y, ,y yy y1 1i iy y计算 的递推公式。 算法：要求给出用已知信息 12专业课件欧拉（Euler）方法是解初值问题的最简单的数值方法。00y)y(xy)f(x,y) )y y, ,( (x x0 00 0y)(x,(x)yy)f(x,9.2 简单的数值方法与基本概念的解y=y(x)代表通过点 的一条称之为微分方程的积分曲线。积分曲线上每一点 的切线的斜率 等于函数 在这点的函数值。 9.2.1 Euler公式初值问题13专业课件Euler法的求解过程: 从初始点P0( 即点(x0,y0)出发，作积分曲线y=y(
6、x)在P0点上切线 ，其斜率为1 10 0P PP P)y,f(x)(xy000y=y(x)x0 xix1yx2P1(x1 , y1)P0Pnxi+1xnP2(x2 , y2)Pi(xi , yi) Pi+1(xi+1 , yi+1)y(x1)y(x2)y(xi)y(xi+1)y(xn)y(x0)14专业课件)x)(xy,f(xyy00001 1x xx x )x)(xy,f(xyy010001这样就获得了P1点的坐标： (x1, y1) 。将y1作为y(x1)的近似值(想象(x1, y1) 在积分曲线y=y(x)上)当 时，得 过点P1(x1,y1)，作积分曲线y=y(x)的切线交直线x=x
7、2于P2点。注意切线 的斜率(近似)为2 21 1P PP P)y,f(x)(xy111直线 方程为:)x)(xy,f(xyy1111)x)(xy,f(xyy121112当 时,得 2xx 由此获得了P2的坐标。直线 的方程为： 1 10 0P PP P2 21 1P PP P15专业课件当 时,得重复以上过程, 对已求得点 ,以 为（近似）斜率作直线 )y,(xPiii)y,f(x)(xyiii)x)(xy,f(xyyiiii1ixx)x)(xy,f(xyyi1iiii1i这样,从x0逐个算出对应的数值解 n21x,x,xn21y,y,y从图形上看,就获得了一条近似于曲线y=y(x)的折线
8、。n321PPPP就获得了一系列的点: P1, P1,Pn。16专业课件y=y(x)x0 xix1yx2y1P0Pnxi+1xny2yiyi+1y(x1)y(x2)y(xi)y(xi+1)y(xn)y(x0)yn微分方程(9.1)的精确解y=y(x)的近似解为： y1， y2， , yn17专业课件注：还可用数值积分法和泰勒展开法推导 Euler公式（略）。Euler公式Euler法的计算公式可以表达为： n0,1,.,i,)y(xy)y,hf(xyy00iii1i ( 9.2 ) 其中， 为常数，i=0,1,ni1ix x h18专业课件nnnnnnn1ny0.2x1.1y)y2xh(yyy
9、解：取h=0.1，根据Euler公式，得1x0 ,y2xyy 1y(0) 例9.1：利用Euler公式求解微分方程的初值问题初值问题有解：2x1y由x0 = 0, y0 =1，代入以上公式，得 y1 = 1.1*y0 - 0.2* x0/y0 = 1.119专业课件课堂练习：计算出 x2 , y2；x3 , y3nnn1ny0.2x1.1yyx0 = 0, y0 =1x1 = 0.1 , y1= 1.120专业课件xnyny(x n)0.11.10001.09540.21.19181.18320.31.27741.26490.41.35821.34160.51.43511.41420.61.5
10、0901.48320.71.58031.54920.81.64981.61250.91.71781.67331.01.78481.7321计算结果比较：初值问题有解：2x1y可以由此公式计算出准确解：y(xn)欧拉法准确值nnn1ny0.2x1.1yy21专业课件y=y(x)的近似解010.1 0.20.30.40.50.60.70.80.9Home11.5222专业课件梯形公式23专业课件9.2.2 梯形公式y)f(x,y x,x1ii1iixxi1iy(x)dxfx,)y(x)y(x1iixxy(x)dxfx,(9.4 ) 改用梯形方法计算其积分项，即 )y(x,f(x)y(x,f(x2x
11、x1i1iiii1i为了提高精度,对方程 的两端在区间 上积分得，24专业课件(9.5)式的右端含有未知的yi+1,它是一个关于yi+1的函数方程,这类数值方法称为隐式方法。相反地,欧拉法是显式方法。 代入(7.4)式,并用近似代替式中即可得到梯形公式 )y,f(x)y,f(x2hyy1i1iiii1i( 9.5 ) 由于数值积分的梯形公式比矩形公式的精度高，所以梯形公式（9.5）比欧拉公式( 9.2 )的精度高。 求解困难Home25专业课件两步Euler公式26专业课件对方程 两端在区间 上积分得 y y) )f f( (x x, ,y y x x, , x x1 1i i1 1i i1i
12、1ixx1i1iy(x)dxfx,)y(x)y(x (9.6 ) 改用中矩形公式计算其积分项，即 )y(x,fxxxy(x)dxfx,ii1i1ixx1i1i代入上式,并用yi近似代替式中y(xi)即可得到)y,2hf(xyyii1i1i( 9.7 ) 9.2.3 两步欧拉公式两步欧拉公式 2个区间27专业课件【注】欧拉方法和梯形方法，都是单步法,其特点是在计算yi+1时只用到前一步的信息yi;而两步欧拉公式 (9.7)中除了yi外,还用到更前一步的信息yi-1,即调用了前两步的信息。 Home28专业课件欧拉法的局部截断误差29专业课件9.2.4 欧拉法的局部截断误差) )y y( (x x
13、y yi ii i1i1iiy)y(xR定义9.1 在yi准确的前提下, 即 时, 用数值方法计算yi+1的误差： 衡量求解公式好坏的一个主要标准是求解公式的精度, 因此引入局部截断误差和阶数的概念。称为该数值方法计算时yi+1的局部截断误差。30专业课件()y2!hy)y(x21i1i (b)-(a），得 )y(xyii)(xyh)y(x)y(x,hf(x)y(xyiiiii1i在欧拉公式中，假定 (近似地相等)，则有(a)而将真解y(x)在xi处按二阶泰勒展开，得 )x,(x()y2!h)(xyh)y(x)y(x1ii2ii1i (b)欧拉公式的截断误差推导：31专业课件定义9.2 若数值
14、方法的局部截断误差为 ，则称这种数值方法的精度阶数是P。)O(h1p 步长(h N 结束。 9.2.6 改进欧拉法算法实现Nh,x,x1057专业课件（2）改进欧拉法的流程图）改进欧拉法的流程图 开开 始始 输输 入入x0, y0,h , N 1 n x0 + h x1 y0+ h f ( x0,y0 ) yp y0+ h f ( x1,yp) yc ( yp+ yc) / 2 y1 输输 出出x1, y1 n + 1 n n = N ? x1 x0 y1 y0 结结 束束 N Y (3) 程序实现(改进欧拉法计算常微分方程初值问题 ) 58专业课件clearx=0, yn=1 %初始化for
15、 n = 1:10yp = yn + 0.1*(yn-2*x/yn); %预测x = x + 0.1;yc = yn + 0.1*(yp-2*x/yp) ; %校正yn = (yp+yc)/2 %平均（再矫正）end本题的精确解为 2x1y(x)改进欧拉法的Matlab程序59专业课件9.3.5 四阶龙格库塔法算法实现Nh,x,x100101yy,xx(1) 计算步骤 输入 使用龙格库塔公式（9.20）计算出y1 输出 ，并使 转到 直至n N 结束。 11y,x60专业课件(2) 四阶龙格库塔算法流程图 开 始 输 入 x0, y0,h , N 1 n x0 + h x1 f(x0,y0 ) k1, f(x0+ h /2 ,y0 + h k1/2 ) k2 f(x0+ h /2 ,y0 + h k2/2 ) k3, f(x1,y0+ h k3) k4 y0+ h (k1+ 2 k2+ + 2 k3+ k4)/6 y1 输 出 x1, y1 n + 1 n n = N ? x1 x0 y1 y0 结 束 n y (3) 程序实现(4阶龙格-库塔法计算常微分方程初值问题) 61专业课件Home62专业课件
