第四章 常微分方程的数值解法 Numerical Solutions to Ordinary Differential Equations
第四章 常微分方程的数值解法 Numerical Solutions to Ordinary Differential Equations
概述 Chapter 4 Initial -value problems for ODE 一阶常微分方程初值问题: d Problem I: dx f(x,y)求函数y=(x)满足: y(x=f(x, y(x) y f(xy)在[a,b]上连续且满足 Lipschitz条件: 彐L,Vy1y2,S.t.|f(x,y1)f(xy2)≤Lyy2 则初值问题 Problem有唯一解 实际工程技术、生产、科研上会出现大量的微分方程问题 很难得到其解析解,有的甚至无法用解析表达式来表示, 因此只能依赖于数值方法去获得徽分方程的数值解。 HUST
概述 ( ) ( ) P roblem Ι : 0 0 d y = f x,y d x y x =y L, y ,y , |f(x, 12 1 2 1 2 ∃∀ ≤ s.t. y )-f(x,y )| L|y | -y !一阶常微分方程初值问题: ! f(x,y)在[a,b]上连续,且满足Lipschitz条件: !实际工程技术 生产 科研上会出现大量的微分方程问题 很难得到其解析解 有的甚至无法用解析表达式来表示 因此只能依赖于数值方法去获得微分方程的数值解 求函数y=y(x)满足 yí(x) =f(x, y(x)) 则初值问题Problem I有唯一解
Chapter 4 Initial -value problems for ODE 微分方程的数值解法: 口不求y=y(X)的精确表达式而求离散点X0x1…X1处的函数值 口设 Problem的解y×)的存在区间是[ab],初始点x=a,取 [a,b]内的一系列节点Xox1…xn。a=x求yory1r…yn HUST
! 微分方程的数值解法 "不求y=y(x)的精确表达式,而求离散点x0,x1,…xn处的函数值 "设Problem I的解y(x)的存在区间是[a,b] 初始点x0=a 取 [a,b]内的一系列节点x0, x1,...,xn a= x0求y0, y1,...,yn
Chapter 4 Initial -value problems for ODE 方法:采用步进式和递推法 将[ab]n等分,a=×1m-1n=m1=)…->Vn=>)n+1- 怎样建立递推公式? Tayo公式 >数值积分法 HUST
! 计算过程 0 0 n+1 n, n n-1 n-2 n-m y = y(x ) y = g(h,x y ,y ,y ,...y ) !!方法 采用步进式和递推法 将[a,b]n等分 a= x0< x1<…< xn =b 步长h= xk=a+kh b-an 0 1 2 nm nm n 1 n 1 y y y yy y y → → →⋅⋅⋅→ → →⋅⋅⋅→ → →⋅⋅⋅ − − + + ! 怎样建立递推公式 !Taylor公式 !数值积分法
4.1欧拉公式 Chapter 4 Initial -value problems for ODE y 狄思想:用向前差商近似代替微商 y(Xo=yo y(x)≈yXn+)y(x) h=Xn+iXn f(xny×n) yxn+,) yX1=yx)+hf(x) yn1=yn+hf(xnyn)—欧拉公式( Euler scheme) Jx%)= n=01,2 Yn+1=y, +hf(xn,y, HUST
!思想:用向前差商近似代替微商. ≈ n+1 n n n+1 n y(x )-y(x ) y'(x ) h=x -x h ⋅⋅⋅⋅⋅⋅ 0 0 n+1 n n n y(x )=y n=0,1,2, y =y +hf(x ,y ) ——欧拉公式 Euler Scheme 4.1 欧拉公式 ∴ ≈ n+1 n n n y(x )- y(x y(x ) f( , x) . ) h ( ) ( ) ' 0 0 y = f x,y y x =y ∴ ≈ n+1 n n n y(x ) y(x )+hf(x ,y(x )) ∴ n+1=y +hf(x ,y n n n y )
几何意义 Chapter 4 Initial -value problems for ODE 1.y(×)过点Pxy0)且 f(x,y 在任意点(X,y)的切线 (o=yo 斜率为f(Xy), 2y(X)在点P0X0y0)的切 线方程为: y=y0+(x82y0)(XX 在切线上取点P1(X1y1 y=yo +hf(xo yo) y1正是Euer公式所求。 3.类似2,过P以f(X1y1为斜率作y(X)的切线,在其上取点 P2(X2y2),依此类推. 4折线P0P1P2Pn…作为曲线y刈)的近似 欧拉折线法 HUST
几何意义 1. y(x)过点P0(x0,y0)且 在任意点(x,y)的切线 斜率为f(x,y), 2. y(x)在点P0(x0,y0)的切 线方程为 3. 类似2 过P1以f (x1,y1)为斜率作y(x)的切线 在其上取点 P2(x2,y2) 依此类推… 4.折线P0 P1 P2 ÖPnÖ作为曲线y(x)的近似 ——欧拉折线法 x p0 p1 p2 p3 p4 x0 x1 x2 x3 x4 y y(x) ( ) ( ) ' 0 0 y = f x,y y x =y 0 00 0 y=y +f(x ,y )(x-x ) 在切线上取点P1 (x1,y1) 1=y + f(x 0 0 0 y h ,y ) y1正是Euler 公式所求
欧拉法(续) Chapter 4 Initial -value problems for ODE 用向后差商近似代替微商 f(,y y(Xn+1)≈yXnX1+1 yxnn1-yX n+1Xn (X0)=Y f(×n+1,y(Xn+1) yx here xn=y(x)+hf(X,+D) y+1=y,+hf(xn+1, n41 隐式欧拉公式 MMx)=y 注:用隐式欧拉法,每一步都需 Yu 1=y, thy nQ12-…解方程(或先解出yn1的显式表达 式),但其稳定性好, 狄用中心差商近似代替微商: y(xn)≈yxn1X yx1)) →f(xnyX) y(xn1)x×n1) 2h 2h Mxo)=y 二步欧拉法 ≌n=m+2fx) 注:计算时,先用欧拉法求出y1,以后再用二步欧拉法计算
!用向后差商近似代替微商 ≈ n+1 n n+1 n n+1 n+1 n y(x )-y(x ) y'(x ) y[x ,x ]= h=x -x h ⋅⋅⋅⋅⋅⋅ 0 0 n+1 n n+1 n+1 y(x )=y n=0,1,2, y =y +hf(x ,y ) ——隐式欧拉公式 欧拉法 续 ∴ ≈ n+ n+ 1 1 n n+1 y(x )-y(x ) f(x , ) y(x ) h ( ) ( ) ' 0 0 y = f x,y y x =y ⇒ n+1 ≈ n n +hf(x , n+1 +1 y(x ) y(x ) y( x ) ) ∴ n 1 + =y +hf( n n+1 n+1 y x , y ) here 注 用隐式欧拉法 每一步都需 解方程 或先解出 yn+1的显式表达 式 但其稳定性好 !用中心差商近似代替微商 ≈ n-1 n+ n+1 n-1 n 1 y(x )-y(x ) y'(x ) y[ , ]= x x 2h ⇒ ≈ n+1 n-1 n n y(x )-y(x ) f(x , ) y(x ) 2h ⋅⋅⋅ n-1 0 0 n+1 n n y(x )=y n=0,1,2, y = +2hf(x ) y , y —— 二步欧拉法 注 计算时 先用欧拉法求出 y1 以后再用二步欧拉法计算
欧拉法(续) Chapter 4 Initial -value problems for ODE 公式 单步否显式否截断误差y(xn+1)yn1 Yn+1=y, +hf(x,y,) 单步显式 yn+1=y, +hf(xn+1,y+1) 单步隐式 yn+1=yn1+2hf(xnyn)二步显式 HUST
欧拉法 续 n+1 n n n y =y +hf(x ,y ) 公式 单步否 显式否 单步 显式 n+1 n n+1 n+1 y =y +hf(x ,y ) 单步 隐式 n+1 n-1 n n y =y +2hf(x ,y ) 二步 显式 截断误差y(xn+1)-yn+1
Chapter 4 Initial 局部截断误差 -value problems for ODE 定义1假设yn=y(x),即第n步计算是精确的前提下称 Rn+1=y(Xn+-)-ym+1为欧拉法的局部截断误差 注:无y=y(x)前提下,称Rn+1为整体截断误差 定义2若某算法的局部截断误差为O(hP+1),称该算法有p阶精度 定理欧拉法的精度是一阶 分析:证明其局部截断误差为O(h2),可通过 Taylor.展开式分析 证明:Euer公式为yn+1=yn+hf(Xnyn) 令yn=y(Xn),下证:y(Xn+)yh+1=O(h2) yn+1=yn+hf(xnyn)=y(Xn)+hf(Xny(Xn)由y(X)=f(x,y(×) x)+h n y(n) Taylo y(xn+1=y(xn +h)===y(x,)+hy(x,)+ y(5) h25∈( 2! n+1 y(xn+D-y h2=O(h2) HUST
局部截断误差 定义1 假设yn=y(xn) ,即第n步计算是精确的前提下 步计算是精确的前提下,称 Rn+1=y(xn+1)-yn+1为欧拉法的局部截断误差 欧拉法的局部截断误差. 定理 欧拉法的精度是一阶 定义2 若某算法的局部截断误差为O(hp+1),称该算法有p阶精度. 注 无yn=y(xn) 前提下 称Rn+1为整体截断误差 分析 证明其局部截断误差为O(h2),可通过Taylor展开式分析 证明 Euler 公式为 n+1=y +hf(x ,y n n n y ) 令yn=y(xn) 下证 y(xn+1)-yn+1 = O(h2) ∵ n+1= +hf(x n n n y y ,y ) = +hf( n n n y(x ) y x , (x )) 由 yí(x) =f(x, y(x)) ' = +h n n y(x ) y (x ) Taylor ' ' n+1 n n n n+1 2 n ' y(x )=y(x +h) === y(x )+hy (x )+ , (x , x y ) 2 ) ( h ξ ξ ∈ ! ∵ '' 2 n+1 ( ) 2 n+1 2 y(x )-y h ) = h =O( y ξ ∴
Chapter 4 Initial 二步欧拉法的局部截断误差 -value problems for ODE y n+1 yn-1 +2hf(X,y,) 定义3假设y=y(×n),Yn=1=yXn-1,称R+1=y(Xn+)-yh+1为 二步欧拉法的局部截断误差 定理隐式欧拉法的精度是一阶,二步欧拉法的精度是二阶 证明:对二步欧拉法进行证明,考虑其局部截断误差, Ayn=y(Xn),y-1=y(X,-1), Yu+1=y 1 +2f( y)=y(x- 1+2hf(x, y(x D)=y(x-1)+2hy(X y(xn+=y(xn +h)===y(x)+hy(xn)+y(Xn)- h,5∈(xn y(xn p=y(x-h)=(x )-hy(xn) y(X)+2m(-h)', nE(Xn-1,X) 3 将上两式左右两端同时相减: y(x)y(x, =2hy(X,)+ y(5)+y(m)13 ∴y(xn+1)-yn+=O(h) 二步欧拉法的局部截断误差为O(h),其精度是二阶。 HUST
二步欧拉法的局部截断误差 二步欧拉法的局部截断误差 定义3 假设yn=y(xn) , yn 1=y(xn 1),称Rn+1=y(xn+1)-yn+1为 二步欧拉法的局部截断误差 二步欧拉法的局部截断误差. n+1= +2 n-1 n n y y y hf(x , ) 定理 隐式欧拉法的精度是一阶 二步欧拉法的精度是二阶 证明 对二步欧拉法进行证明 考虑其局部截断误差 令yn=y(xn) , yn 1=y(xn 1), = +2hf n n -1 n y(x ) y (x , (x )) ' = + n-1 n y( x ) y (x ) 2 h Taylor ' 2 '' n+1 n n n n n+ ' 1 ' ' 3 n y(x )=y(x +h) === y(x )+hy (x ) h y (x ) 2! + + ,( y ( ) ) h x ,x ξ ξ ∈ 3! ∵ n+1 n+ 3 1 ∴ y(x )-y =O(h ) n+1= +2 n-1 n n y y y hf(x , ) ' n-1 n n n 2 n-1 '' 3 n n ' '' (-h) y (x ) 2 y ( y(x )=y( ) x -h)=y(x )-hy (x )+ + , (x , x ) (-h) ! η η ∈ 3! 将上两式左右两端同时相减 将上两式左右两端同时相减 ''' ''' ' n-1 () () 3 n 1+ n 3 y(x )-y(x ) 2hy (x = )+ h y y ξ η + ∴ ! 二步欧拉法的局部截断误差为 二步欧拉法的局部截断误差为O(h3) 其精度是二阶
