1 大学数学实验 Experiments in Mathematics 实验 4 常微分方程数值解 清华大学数学科学系 为什么要学习微分方程数值解 ? 微分方程是研究函数变化规律的重要工具，有着广泛 的应用。如 : 物体的运动 , 电路的电压 , 人口增长的预测 ? 许多微分方程没有解析解，数值解法是求解的重要手 段，如 2 , dy yx dx = + () () xt axy yt axy by =? =? & & 实验4的基本内容 3. 实际问题用微分方程建模，并求解 2. 龙格 -库塔方法的 MATLAB实现 *4. 数值算法的收敛性、稳定性与刚性方程 1. 两个最常用的数值算法 : ? 欧拉（ Euler）方法 ? 龙格 -库塔（ Runge-Kutta）方法 实例1 海上缉私 海防某部缉私艇上的雷达发现正东方向 c海里处有一艘走私船 正以速度 a向正北方向行驶，缉私艇立即以最大速度 b(>a)前往 拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终 指向走私船。 ? 建立任意时刻缉私艇位置及 航线的数学模型 ,并求解 ; ? 求出缉私艇追上走私船的时间。 a 北 b c 艇船 实例1 海上缉私 建立坐标系如图 : t=0 艇在 (0, 0), 船在 (c, 0); 船速 a, 艇速 b 时刻 t 艇位于 P(x, y), 船到达 Q(c, at) 模型 : 0 y x c R(c,y ) α Q(c,at) P(x,y) b cos , sin dx dy bb dt dt α α== 22 22 () ()( ) () ()( ) dx b c x dt cx aty dy b at y dt cx aty ?? = ? ?+? ? ? ? ? = ? ?+? ? 由方程无法得到 x(t), y(t)的解析解 需要用数值解法求解 “常微分方程初值问题数值解”的提法 00 '(,),()yfxyyx y= =设 的解y=y(x)存在且唯一 不求解析解 y = y(x)(无解析解或求解困难 ) 12 n xx x< <<<LL而在一系列离散点 () ( 1,2,) n yx n= L n 求 的近似值y 通常取等步长 h 0n x xnh= + ? 0 x 0 y )(xyy = y x y 1 y 2 y n )( n xy 1 x 2 x n x 2 y P 0 x 0 x 1 x 2 x 3 x y=y(x) y 0 P 1 P 2 P 3 欧拉方法 00 '(,),()y fxy yx y== 基本 思路 1 [, ] nn xx + 在小区间 上 1 (, ) [ , ] nn fxy xx + 中的x取 内的某一点 1 '[( ) ()]/, nn yyx yxh + = ? 11 ( ) ( ) ( , ( )), [ , ] nn n yx yx hfxyx x x x ++ =+ ∈ x取不同点 各种欧拉公式 n x x取左端点 1 ()() (,() nnnn yx yx hfx yx + =+ 1 (, ), 0,1, nn nn yyhfxyn + =+ =L 向前欧拉公式 显式公式 11 (), ( ) nnn n yyxy yx ++ ≈≈ P 3 欧拉方法 11 ( ) ( ) ( , ( )), [ , ] nn n yx yx hfxyx x x x + + = +∈ 11 (), ( ) nnn n x yyxy yx + + ≈≈取右端点， 111 (, ),0,1, nn nn yyhfxyn ++ = +=L 向后欧拉公式 隐式公式 y P 0 x 0 x 1 x 2 x 3 x y 0 y=y(x) P 1 P 2 n+1 y右端 未知，需迭代求解 (0) 1 (, ) nn nn yyhfxy + =+初值 (1) () 111 (, 0,1,2, , 0,1,2, kk nn nn yyhfxy kn + +++ =+ ==LL () 11 lim k nn k yy ++ →∞ = 向前欧拉公式 向后欧拉公式 1 (, ) nn nn yyhfxy + =+ 111 (, ) nn nn yyhfxy + ++ =+ 二者平均得到 梯形公式 111 [( , ) ( , )], 0,1, 2 nn nn nn h yy fxyfxyn ++ =+ + = L 仍为隐式公式，需迭代求解 将梯形公式的迭代过程简化为两步 1 (, ) nn nn yyhfxy + =+ 预测 111 [( , ) ( , )], 0,1, 2 nn nn nn h yy fxyfxyn +++ =+ + = L 校正 改进欧拉公式 龙格-库塔方法 11 ( ) ( ) ( , ( )), [ , ] nn n yx yx hfxyx x x x + + =+ ∈ ? 向前，向后欧拉公式： 用 [x n , x n+1 ]内１个点的导 数代替 f(x, y(x)) ? 梯形公式，改进欧拉公式： 用 [x n , x n+1 ]内２个点导数的平均值代替 f(x, y(x)) 龙格-库塔方法的基本思想 在 [x n , x n+1 ]内多取几个点，将它们的导数加权平 均代替 f(x, y(x))， 设法构造出 精度更高 的计算公式。 常用的 (经典 ) 龙格 —库塔 公式 不足 :收敛速度较慢 1 1 1 2221 1 1 (, ) (, ) (, ),3,4, L nn i i nn nn i ininiij j yyhk kfxy kfxchychk kfxchychak i L λ + = ? = ? =+ ? ? ? = ? =+ + ? ? ??? ? ? =+ + = ? ? ∑ ∑ L L级 ?阶龙格 -库塔 方法的 一般形式 11234 1 21 32 43 (22 ) 6 (, ) (/2, /2) (/2, /2) (, ) nn nn nn nn nn h yy kkkk kfxy kfxhyhk kfxhyhk kfxhyhk + ? =+ + + + ? ? = ? ? =+ + ? ? =+ + ? =++? ? ? 1,10,1,, 1 11 =≤≤= ∑∑ ? == i j iji L i iijii acac λλ 满足 使精度尽量高 4级 4阶 微分方程组和高阶方程初值问题的数值解 欧拉方法和龙格 -库塔方法可直接推广到微分方程组 0000 (,,) (,,) () ,() yfxyz zgxyz yx y zx z ′= ? ? ′= ? ? = = ? 向前欧拉公式 ),,( yyxfy ′=′′ 1 yy= ? ? ? =′ =′ ),,( 212 21 yyxfy yy 高阶方程需要先降阶为一阶微分方程组 ? ? ? ? ? = += += + + L,2,1,0 ),,( ),,( 1 1 n zyxhgzz zyxhfyy nnnnn nnnnn 3 龙格 —库塔方法的 MATLAB 实现 T n T n fffxxxxtxxtftx ),(,),(,)(),,()( 1100 LL& ==== [t,x]=ode23(@f, ts, x0) 3级 2阶龙格 -库塔公式 [t,x]=ode45(@f, ts, x0) 5级 4阶龙格 -库塔公式 f是待解方程写成 的函数 m文件： function dx=f(t,x) dx=[f1; f2;… ; fn]; ts = [t0,t1, …,tf] 输出指定时刻 t0,t1, …,tf 的函数值 ts = t0:k:tf 输出 [t0,tf] 内等分点处的函数值 x0为函数初值 (n维 ) 输出 t=ts, x为相应函数值 (n维 ) 缺省精度(相对误差10 -3 ，绝对误差10 -6 ), 计算步长按精度要求自动调整. 实例 1 海上缉私(续) 模型的数值解 22 22 () ()( ) () ()( ) dx b c x dt cx aty dy b at y dt cx aty ?? = ? ?+? ? ? ? ? = ? ?+? ? (0) 0, (0) 0xy= = 0 y x c (x(t),y(t)) a b 设 : 船速 a=20 (海里 /小时 ) 艇速 b=40 (海里 /小时 ) 距离 c=15 (海里 ) 求 : 缉私艇的位置 x(t),y(t) 缉私艇的航线 y(x) x0=[0 0]; a=20;b=40;c=15; [t,x]=ode45(@jisi,ts,x0,[],a,b,c); %exact solution x1=c y1=a*t; %output t,x(t),y(t) [t,x,y1] %draw x(t),y(t) plot(t,x),grid, gtext('x(t)','FontSize',16), gtext('y(t)','FontSize',16),pause %draw y(x): the position of tatch js plot(x(:,1),x(:,2),'r*'),grid xlabel('x','FontSize',16), ylabel('y','FontSize',16) %Creat the function for jisi.m %Let x(1)=x, x(2)=y function dx=jisi(t,x,a,b,c) s=sqrt((c-x(1))^2+(a*t-x(2))^2); dx=[b*(c-x(1))/s;b*(a*t-x(2))/s]; MATLAB 6.5.1.lnk jisi.m, seajisi.m 实例 1 海上缉私(续) 模型的数值解 a=20, b=40, c=15 走私船的位置 x 1 (t)= c=15 y 1 (t)=at=20t t=0.5时缉私艇追上走私船 0 5 10 15 20 0 2 4 6 8 10 x y 缉私艇的航线 y(x) 9.997915.00460.50 8.027314.74510.45 6.177313.98060.40 4.555212.81700.35 3.200511.34960.30 2.11789.67050.25 1.28997.85150.20 0.69065.94450.15 0.29243.98540.10 0.06981.99840.05 000 y(t)x(t)t 10.0 9.0 8.0 7.0 6.0 5.0 4.0 3.0 2.0 1.0 0 y 1 (t) 实例 1 海上缉私(续) 模型的数值解 设 b， c不变， a变大 为 30, 35, … 接近 40, 观察解的变化 : a=35, b=40, c=15 t=? 缉私艇追上走私船 56.055.948614.98661.6 52.552.014615.00231.5 49.048.018315.01171.4 45.544.016514.99961.3 42.040.016414.99861.2 ………… 17.512.083013.75510.5 14.08.275512.53840.4 10.54.828310.52400.3 7.02.13087.59280.2 3.50.50583.95610.1 0000 y 1 (t)y(t)x(t)t 累 积误差较大 提高精度! 实例 1 海上缉私(续) 模型的数值解 a=35, b=40, c=15 opt=odeset('RelTol',1e-6, 'AbsTol',1e-9); [t,x]=ode45(@jisi,ts,x0,opt); t=1.6时缉私艇追上走私船 0 5 10 15 20 0 20 40 60 x y 缉私艇的航线 y(x) 判断 “追上 ”的有效方法 ? 56.055.99993115.0000201.6 52.552.00000514.9999981.5 49.048.00000514.9999931.4 45.544.00000514.9999631.3 42.040.00000514.9996161.2 ………… 17.512.07534413.7539740.5 14.08.26984012.5394540.4 10.54.82930810.5219210.3 7.02.1306787.5928220.2 3.50.5058133.9561040.1 0000 y 1 (t)y(t)x(t)t 4 实例 1 海上缉私(续) 模型的解析解 2222 )()( )( , )()( )( yatxc yatb dt dy yatxc xcb dt dx ?+? ? = ?+? ? = )( xc yat dx dy ? ? = aty dx dy xc =+? )( 2 2 () dy dt cx a dx dx ?= ds b dt = 22 () ()ds dx dy=+ 2 2 2 () 1() dy a dy cx dx b dx ?=+ dx dy p =令 b a kpk dx dp xc =+=? ,1)( 2 (0) 0p = 实例 1 海上缉私(续) 模型的解析解 2 1)( pk dx dp xc +=? (0) 0p = 2 1() k cx pp c ? ? ++= 2 1() k cx pp c ? +?= 1 [( ) ( ) ] 2 kk cx cx p cc ? ?? =? /1kab= < (0) 0y = 11 2 11 [() ()] 21 1 1 kk ccx cx kc y kc kc k +? ?? =?+ + ?? 缉私艇的航线 y(x)的解析解 x=c时 缉私艇追上走私船的 y坐标 缉私艇追上走私船的时间 : 1 22 bc t ba = ? a=20, b=40, c=15 → t 1 =0.5 a=35, b=40, c=15 → t 1 =1.6 222 1 kc abc y kba == ? ? 微分方程数值解法的误差分析 数值解法 : 计算微分方程精确解 y(x n )的近似值 y n 按照步长 h一步步计算 , 每步都有误差 ; 每一步的误差会逐步积累, 称 累积误差 . 讨论计算一步出现的误差 假定 y n = y(x n ) , 估计 y n+1 的误差 : y(x n+1 )- y n+1 局部截断误差 误差分析 估计欧拉公式的局部截断误差 y(x n+1 )在 x n 处作 Taylor展开 : )()( 2 )()()( 3 2 1 hOxy h xyhxyxy nnnn +′′+′+= + 向前欧拉公式 1 (, ) nn nn yyhfxy + = + y n = y(x n ) )()())(,()( 1 nnnnnn xyhxyxyxhfxyy ′+=+= + 2 32 111 () ()() () 2 nnn n h T yx y y x Oh Oh +++ ′′=?= += 局部截断误差主项为 )( 2 2 n xy h ′′ 误差分析估计欧拉公式的局部截断误差 )()( 2 )()()( 3 2 1 hOxy h xyhxyxy nnnn +′′+′+= + 向后欧拉公式 ),( 111 +++ += nnnn yxhfyy 2 32 111 () ()() () 2 nnn n h Tyx y yxOhOh +++ ′′=?=? += 局部截断误差主项为 )( 2 2 n xy h ′′? 梯形 公式 向前、向后欧拉公式的平均 3 43 1 () () () 12 nn h TyxOhOh + ′′′=? + = x n x n+1 x y P n A Q B 误 差 分 析 算法 精度的阶 的定义 一个算法的局部截断误差为 O(h p+1 ) 该算法具有 p阶精度 局部截断误差 精度 向前欧拉公式 O(h 2 ) 1阶 向后欧拉公式 O(h 2 ) 1阶 梯形公式 O(h 3 ) 2阶 改进欧拉公式 O(h 3 ) 2阶 经典龙格-库塔公式 O(h 5 ) 4阶 5 实例 2 弱肉强食 问 题 自然界中同一环境下两个种群之间的生存方式 相互竞争 相互依存 弱肉强食 弱肉 强食 种群甲靠丰富的自然资源生存 食饵(Prey) 种群乙靠捕食种群甲为生 捕食者(Predator) 两个种群的数量如何演变? 实例 2 弱肉强食 模型 食饵 (甲 ) 的密度 x(t), 捕食者 (乙 )的密度 y(t) 0,/ >= rrxx& 甲独立生存的增长率 r 0,/ >?= aayrxx& 乙使甲的增长率减小, 减小量与 y 成正比 0,/ >?= ddyy& 乙独立生存的死亡率 d 0),(/ >??= bbxdyy& 甲使乙的死亡率减小, 减小量与 x成正比 00 () () (0) , (0) xrayxrxaxy ydbxydybxy xxyy =? =? ? ? =? ? =? + ? ? == ? & & Volterra模型 x(t), y(t)无解析解 实例 2 弱肉强食 模型的数值解 ? ? ? ? ? == +?= ?= 00 )0(,)0( yyxx bxydyy axyrxx & & 00 1, 0.5 0.1, 0.02 25, 2 rd ab xy == == == 0 5 10 15 0 20 40 60 80 100 x(t) y(t) 0 20 40 60 80 100 0 5 10 15 20 25 30 x y 猜测 x(t), y(t)是周期函数 ; y(x)是封闭曲线 数值积分计算一个周期的平均值： 25, 10xy≈≈ MATLAB 6.5.1.lnk shier.m, shier1.m ybxd xayr dy dx )( )( ?? ? = 实例 2 弱肉强食 模型的解析解 ? ? ? ??= ?= ybxdy xayrx )( )( & & ceyex ayrxxd = ?? ))((相轨线 dy y ayr dx x bxd ? = +? c由初始条件确定 相轨线是封闭曲线 (c在一定范围内 ) 求 x(t),y(t)一周期的平均值 : yx , 可以 证明 ybxdty )()( ??=& bdyytx /)/()( += & ∫ = T dttx T x 0 )( 1 b d x = x(t), y(t)是周期函数 (周期记作 T) ) )0(ln)(ln ( 1 b dT b yTy T + ? = 实例 2 弱肉强食 模型的解析解 r y a =x(t),y(t)一周期的平均值 : b d x = 02.0,1.0,5.0,1 ==== badr 10,25 == yx r ~食饵增长率 a ~捕食者对食饵的捕获能力 d ~捕食者死亡率 b ~食饵对捕食者的喂养能力 结果解释 ? ? ? ??= ?= ybxdy xayrx )( )( & & 与计算结果同 ↑?↓↑ yar , ↑?↓↑ xbd , 既相互制约 又相互依存 欧拉和龙格 -库塔方法的共同点 :只用 y n 计算 y n+1 1 00 (,,) () nn nn yyhxyh yyx ? + =+? ? = ? 单步法 步长 h→0时数值解 y n 无限接近解析解 y(x n ) 收敛性 单步法有 p阶精度 (局部截断误差 O(h p+1 )) (, ,) (, ,) ( 0)xyh xyh Ly y L??? ≤? > 向前欧拉公式 改进欧拉公式 4阶龙格-库塔公式 数值算法的收敛性和稳定性 (显式)单步法: 整体误差 e n =O(h) e n =O(h 2 ) e n =O(h 4 ) 单步法收敛 整体误差 e n =y(x n )-y n =O(h p )→0 (h→0) 6 数值算法的收敛性和稳定性 稳定性 计算中舍入误差不会随步数的增加无限增大 L,2,1, =≤ + k nkn εε y n 的误差 ε n ),( yxfy =′ λ>0→微分方程稳定 ),( 1 nnnn yxhfyy += + 向前欧拉公式 向后欧拉公式 11nn n yyhyλ ++ =? 经典龙格-库塔公式 2.785/h λ≤ 算法稳定 ))(,())(,(),( ******** yyyxfxxyxfyxfy yx ?+?+=′ 0, >?=′ λλyy (特征根 -λ) nn h ελε )1( 1 ?= + nk n ε ε + ≤ 11hλ? ≤ 2/h λ≤ 1 1 1 nn h ε ε λ + = + h任意 x yce λ? = n yh )1( λ?= 刚性现象与刚性方程 bxaxrktfrxxkx ==>=++ )0(,)0(),0,()( &&&& 振动系统或电路系统的数学模型现象 k=2000.5, r=1000, a=1, b=-1999.5, f(t)=1 瞬态解与稳态解 e -2000t ~快瞬态解 e -t/2 ~慢瞬态解 计算到 t=0.005时已衰减到 4.5×10 -5 计算到 t=20时才衰减到 4.5×10 -5 精度达到 10 -4 需算到 t=20(由慢瞬态解 λ=1/2决定 ) 选取步长 h由快瞬态解 λ=2000决定 h<2.785/2000=0.0014 龙格 -库塔公式 t=20需 14286步 快 、 慢瞬态解的特征根相差悬殊 2000 / 2 () 1 tt xt e e ?? = ?+ 刚性现象( Stiff) 求稳 态解 刚性现象与刚性方程 刚性 方程 振动、电路及化学反应中的线性常系数方程组 n RxtfAxtx ∈+= ),()(& A的特征根 λ k (k=1,2,… ,n) 的实部 Re(λ k )<0 max Re( ) min Re( ) k k k k s λ λ = ~快瞬态解大)Re( k λ ~慢瞬态解 小)Re( k λ ~刚性比 快瞬态解 →步长充分小 慢瞬态解 →积分区间长 s>10 刚性非线性方程组 线性常系数方程组 传统的数值方法无能为力 线性化方法 刚性方程 刚性现象与刚性方程 刚性方程的 MATLAB求解 ode23, ode45 解刚性方程的困难 步长自动变小 计算时间很长 求解刚性方程的命令 : ode23s, ode15s 等 (用法相同) 11 2 66 212 12 2 (10 1) (10 2) (0) 10 / 4, (0) 10 / 4 1/ 2 xx x x xx =+? ? =? + ? + ? ? ==? ? & & 例 特征根 λ 1 =-1,λ 2 =-10 6 刚性比 s=10 6 6 6 6 10 1 66 10 2 10 () ( 1) 4 10 (10 1) () ( 1) 42 tt tt xt e e xt e e ?? ?? ? =+? ? ? ? + ? =? + + ? ? 解析解 MATLAB 6.5.1.lnk Stiff.m,stiff1.m 布置实验 目的 1.用 MATLAB软件掌握求微分方程数值解的方 法，并对结果作初步分析； 2.通过实例学习用微分方程模型解决简化的实际问 题。 内容 课上布置，或参见网络学堂