大学酸学实 大学数学实验 为什么要学习微分方程数值解 微分方程是研究函数变化规律的重要工具,有着广泛 的应用。如: Experiments in Mathematics 物体的,电路的电压人口增长的顶 ·许多微分方程没有解析解,数值解法是求解的重要手 实验4常微分方程数值解 段,如 dy --axx 清半大总教科導系 y(n)=axy-by 实验4的基本内容 实例1海上缉私 海防某部缉私艇上的达发现正东方向c海里处有一艘走 1.两个最常用的数值算法: 正以速度a向正北方向行驶,盘私艇立即以最大遠度b>a) 欧拉(Euer)方法 拦截。如果用雷达进行跟踪时,可保持私艇的速度方向始终 龙格库塔( Runge-Kutt)方法 向走私部 2.龙格库塔方法的 MATLAB实现 建立任意时刻缉私艇位量及 航线的数学模型并求解; 3.实际问题用微分方程建模,并求解 求出缉私艇追上走私船的时间 4数值算法的收敛性、稳定性与刚性方程 (学静学实鉴 学学实纷 实例1海上缉私 “常微分方程初值问题数值解”的提法 立系如圆:P0艇在(0,0),船在(c,0);船速a,艇速b 时刻t艇位于Px,y),船到达Qc,a 设y=f(x,y)2y(x)=y的解y=y(x)存在且唯一 不求解析解y=y(x)(无解析解或求解困难) -= bcos -=bsin a 而在一系列高散点<x<…<x dt )的近似值ya(n=1,2 dx c-x)2+(at-y) 通常取等步长h y=y(x) y(x, 由方程无法得到x(,y(的解析解 需要用数值解法求解 No 133
1 大学数学实验 Experiments in Mathematics 实验4 常微分方程数值解 清华大学数学科学系 为什么要学习微分方程数值解 • 微分方程是研究函数变化规律的重要工具，有着广泛 的应用。如: 物体的运动, 电路的电压, 人口增长的预测 • 许多微分方程没有解析解，数值解法是求解的重要手 段，如 2 , dy y x dx = + ( ) ( ) x t axy y t 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 b b dt dt = = α α 2 2 2 2 ( ) ( )( ) ( ) ( )( ) dx b c x dt c x at y dy b at y dt c x at y ⎧ − = ⎪ ⎪ −+− ⎨ − ⎪ = ⎪ − +− ⎩ 由方程无法得到x(t), y(t)的解析解 需要用数值解法求解 “常微分方程初值问题数值解”的提法 0 0 设 的解y=y(x)存在且唯一 y f xy yx y ' ( , ), ( ) = = 不求解析解 y = y(x) (无解析解或求解困难) 1 2 n 而在一系列离散点 xx x < <<< L L ( ) ( 1, 2 , ) n 求 的近似值y yx n n = L 通常取等步长h n 0 x = x nh + • 0 x 0 y y = y(x) y x y1 y2 yn ( ) n y x 1 x 2 x n x
大学静学实扮 欧拉方法y=f(xy)y(x)=y 欧拉方法y(xn)=y(xn)+f(x,y(x)x∈[x2,x 基本在小区间x,x上y=[y(xm)-y(xn)/h x取右端点,yn=y(x)ym≈H(xn)几 思路f(x,y)中的x取x,x内的某一点 yn1=yn+(xn1,yn)=O,1…,向后歇拉公式 隐式公式 M=xr,)+hf(x,Mr)),xElruru-I 右端y。1未知,需迭代求解 x取左端点x,几 各种欧拉公式 Jour) 初值y1=yn+bf(xn,yn) M(rm)=xr,)+f(rn, y(r, ) ymt=y, +hf(r,I,ym)) yn=y(x)y≈(xn) ym=+的(xy)=… y 向前欧拉公式 显式公式 向前欧拉公式 向后欧拉公式 龙格一库塔方法 yu.l=v.+ hf(xu, yu) y(m =y(x)+hf(x,y(x)),xELx xnI 二平站梯形公 劑,向看吹公式 f(r,y)+f(xm,yu)l, n=0,1 xnxa+11个的寻代警几x,y(x) 仍为隐式公式,需选代求解 御形公。欧公 改进欧拉公式将梯形公式的选代过程简化为两步 用[pxal肉2个羸款的平代兽几(x,y(x) n=yn+bf(xn,y)预测 龙格一库塔方法的基本思想 y=y+,f(x,y)+f(xa1,),n=0,1, 在xnxm内多取几个点,将它们的导数加权平 均代替fx,yx),设法构造出精度更高的计算公式 校正 (学静学实鉴 (大学数学实验) 龙格库塔 分方翟组和痛阶方视物值向题的歙 K,=f(r,y) 方法的 一般形式 k2=f(r,+c,h,y,+c,hk,) 欧拉方法和龙格库塔方法可直接推广到微分方程组 y′=f(x,y,=) yn+的f(xn,yn,三n) k=f(xn+ch+c小∑ak,1=3,4… λ,c,a满足∑A,=1,0sc,≤1,∑a=1使精度尽量高 y(x0)=y,(x0)==0 向前歌拉公式 n=y+(k1+2k2+2k+k) 高阶方程需要先降阶为一阶微分方程组 常用的(经 龙格一库塔 Ey= f(x+h/2,y+hk1/2) y”=f(xyy)y y2=f(r,y,y2) 不足收敛遠度较慢{=八(x+,+航
2 y P0 x0 x1 x2 x3 x y=y(x) y0 P1 P2 P3 欧拉方法 0 0 y f xy yx y ' ( , ), ( ) = = 基本 思路 1 [, ] n n x x 在小区间 上+ 1 (, ) [ , ] n n f xy x x 中的x取 内的某一点 + 1 ' [ ( ) ( )]/ , n n y yx yx h = + − 1 1 ( ) ( ) ( , ( )), [ , ] n n nn y x y x hf x y x x x x + + =+ ∈ x取不同点 各种欧拉公式 n x取左端点x 1 ( ) ( ) ( , ( )) n n nn y x y x hf x y x + = + 1 ( , ), 0,1, n n nn y y hf x y n + =+ = L 向前欧拉公式 显式公式 1 1 ( ), ( ) n nn n y yx y yx ≈ ≈ + + P3 欧拉方法 1 1 ( ) ( ) ( , ( )), [ , ] n n nn y x y x hf x y x x x x + = + ∈ + 1 1 ( ), ( ) n nn n x y yx y yx 取右端点， ≈ ≈ + + 1 11 ( , ), 0,1, n n nn y y hf x y n + ++ = + = L 向后欧拉公式 隐式公式 y P0 x0 x1 x2 x3 x y0 y=y(x) P1 P2 n+1 右端 未知，需迭代求解 y (0) 1 (, ) n n nn y y hf x y 初 值 + = + ( 1) ( ) 1 11 (, ) 0,1, 2, , 0,1, 2, k k n n nn y y hf x y k n + + ++ = + = = L L ( ) 1 1 lim k n n k y y + + → ∞ = 向前欧拉公式 向后欧拉公式 1 (, ) n n nn y y hf x y + = + 1 11 (, ) n n nn y y hf x y + = + + + 二者平均得到梯形公式 1 11 [ ( , ) ( , )], 0,1, 2 n n nn n n h y y fx y fx y n + ++ =+ + = L 仍为隐式公式，需迭代求解 将梯形公式的迭代过程简化为两步 1 (, ) n n nn y y hf x y + = + 预测 1 11 [ ( , ) ( , )], 0,1, 2 n n nn n n h y y fx y fx y n + ++ =+ + = L 校正 改进欧拉公式 龙格-库塔方法 1 1 ( ) ( ) ( , ( )), [ , ] n n nn y x y x hf x y x x x x + =+ ∈ + • 向前，向后欧拉公式： 用[xn, xn+1]内１个点的导 数代替 f(x, y(x)) • 梯形公式，改进欧拉公式： 用[xn, xn+1]内２个点导数的平均值代替 f(x, y(x)) 龙格-库塔方法的基本思想 在[xn, xn+1]内多取几个点，将它们的导数加权平 均代替 f(x, y(x))，设法构造出精度更高的计算公式。 常用的(经典) 龙格—库塔 公式 不足:收敛速度较慢 1 1 1 2 2 21 1 1 (, ) (, ) ( , ), 3,4, , L n n ii i n n n n i i n i n i ij j j y yh k k fx y k f x c h y c hk k f x ch y ch a k i L + λ = − = ⎧ ⎪ = + ⎪ ⎪ = ⎪ ⎨ =+ + ⎪⋅⋅⋅ ⎪ ⎪ =+ + = ⎪ ⎩ ∑ ∑ L 龙格-库塔 L级?阶 方法的 一般形式 1 1 2 34 1 2 1 3 2 4 3 (22 ) 6 (, ) ( / 2, / 2) ( / 2, / 2) (, ) n n n n n n n n n n h y y k k kk k fx y k f x h y hk k f x h y hk k f x h y hk + ⎧ =+ + + + ⎪ ⎪ = ⎪ ⎪ ⎨ =+ + ⎪ =+ + ⎪ ⎪ =++ ⎪ ⎩ , , 1, 0 1, 1 1 1 1 ∑ = ≤ ≤ ∑ = − = = i j i ij L i i i ij i λ c a 满足 λ c a 使精度尽量高 4级4阶 微分方程组和高阶方程初值问题的数值解 欧拉方法和龙格-库塔方法可直接推广到微分方程组 0 000 (,,) (,,) ( ) ,( ) y f xyz z gxyz yx y zx z ⎧ ′ = ⎪ ⎨ ′ = ⎪ ⎩ = = 向前欧拉公式 y′′ = f (x, y, y′) 1 y y = ⎩ ⎨ ⎧ ′ = ′ = ( , , ) 2 1 2 1 2 y f x y y y y 高阶方程需要先降阶为一阶微分方程组 ⎪ ⎩ ⎪ ⎨ ⎧ = = + = + + + 0,1,2,L ( , , ) ( , , ) 1 1 n z z hg x y z y y hf x y z n n n n n n n n n n
龙格席塔方法的 MATLAB实親 实例1海上私() ()=f(t,x),x(t0)=x0,x=(x1,…xn),f=(f1,…fn)y 模型的数值解 [x=0de23(@t,ts,xo)32阶龙格库塔公式 x=ode45(@t,ts,x)5氟4阶龙格库塔公式 +(ar 是待解方程写成 function dx=f(tx) xin,(o) 的函数m文件 √(c-x)2+(at-y ts=[to,t1,…,t输出指定时刻to,t,….,t的函数值 ts to: k: tf 输出[to切内等分点处的函数值 船速a=20(海里小时) 求:私晨的世置x(2y(0 艇速b=40(海里/小时) 私的航能y(x) X0为函数初值(n维)输出t=ts,x为相应函数值(n维) 缺省精度(相对误差103,绝对误差106, 距高c=15(海里 计算步长按精度要求自动调整 output t,x(D), y(t) ∥寥例1海上私()型的激 a=20b=40,c=15; a=20,b=40,c=15 x =ode 45(@jisi, ts, xo, D, a, b, c) plot(t, x),grid, gtext(,x(t), FontSize, 16 走私船的位量 gtext('y(ty, FontSize, 16),pause x1(D)=c=15 0103.98540.29242.0 %draw y(x): the position of match js %Creat the function for jisi. m plot(x( 1), x( 2), T"), grid 055945410690630 缉私艇的航线v(x) 0207851519940 sqrt(c-x(1)^2+a·-x(2)2 0.301134963.20056.0 dx=|b°(c-x(1)/sb°a·tx(2)s 0351217045570 =0.5时缉私艇追上走私船0451474518.027390 被的值解 学学实 实创1 模型的数值解学醇学实妇 实1 海上() yoo v,a 海上輯私() 设,不变,凌变大03%58 为30,35接近40,03105240482831105 观察解的变化 041253482755140 地的航yx 0412539454826984014 a=35,b=4,c=150513751120830175 0.51375397412.075344175 ?缛私艇追上走私船□121498940144:20 214996100010 13「14996441514551 13|149941001455 积慑戆大 1415011480183490 提高精度 产6时缉私艇追上走私船141900 515003520146525 6149865594865 判“地上的有欺亦海?|Ld1sms9g
3 龙格—库塔方法的 MATLAB 实现 T n T n x(t) f (t, x), x(t ) x , x (x , x ) , f ( f , f ) & = 0 = 0 = 1 L = 1 L [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 海上缉私(续) 模型的数值解 2 2 2 2 ( ) ( )( ) ( ) ( )( ) dx b c x dt c x at y dy b at y dt c x at y ⎧ − = ⎪ ⎪ −+− ⎨ − ⎪ = ⎪ − +− ⎩ x y (0) 0, (0) 0 = = 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 走私船的位置 x1(t)= c=15 y1(t)=at=20t t=0.5时缉私艇追上走私船 0 5 10 15 20 0 2 4 6 8 10 x y 缉私艇的航线y(x) 0.50 15.0046 9.9979 0.45 14.7451 8.0273 0.40 13.9806 6.1773 0.35 12.8170 4.5552 0.30 11.3496 3.2005 0.25 9.6705 2.1178 0.20 7.8515 1.2899 0.15 5.9445 0.6906 0.10 3.9854 0.2924 0.05 1.9984 0.0698 0 0 0 t x(t) y(t) 10.0 9.0 8.0 7.0 6.0 5.0 4.0 3.0 2.0 1.0 0 y1(t) 实例1 海上缉私(续) 模型的数值解 设b，c不变，a变大 为30, 35, …接近40, 观察解的变化 : a=35, b=40, c=15 t=? 缉私艇追上走私船 1.6 14.9866 55.9486 56.0 1.5 15.0023 52.0146 52.5 1.4 15.0117 48.0183 49.0 1.3 14.9996 44.0165 45.5 1.2 14.9986 40.0164 42.0 … … … … 0.5 13.7551 12.0830 17.5 0.4 12.5384 8.2755 14.0 0.3 10.5240 4.8283 10.5 0.2 7.5928 2.1308 7.0 0.1 3.9561 0.5058 3.5 0 0 0 0 y1 t x(t) y(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) 判断“追上”的有效方法? 1.6 15.000020 55.999931 56.0 1.5 14.999998 52.000005 52.5 1.4 14.999993 48.000005 49.0 1.3 14.999963 44.000005 45.5 1.2 14.999616 40.000005 42.0 … … … … 0.5 13.753974 12.075344 17.5 0.4 12.539454 8.269840 14.0 0.3 10.521921 4.829308 10.5 0.2 7.592822 2.130678 7.0 0.1 3.956104 0.505813 3.5 0 0 0 0 y1 t x(t) y(t) (t)
大学静学实扮 寥例1海上私(蟆)模型的解斬解 寥例1海上私()模型的解析解 √e-x)+(a-y3d√e-x)+(a-y p(0)=0 -x)*+y=ar 4(c-x =a业 y(0)=0 舞的航能y(x)的靠折 (c-x =c时 1-k2=b2-a和地上竞船的质 ds=√dx)2+(dy)2 缉私艇追上走私船的时间:1 P a=20b=40,=15→1=05=35,b=40,c=15→1=16 分方翟歙官解油的误分析 误差分析估计欧拉公式的局部截断误差 数值解法:计算微分方程精确解y(xn)的近似值y xm)在xn处作 Taylor展开 按照步长h一步步计算每步都有误差; Mxm)=y(r,)+hy(r, )+y(.)+O(h 一步的误差会逐步积累,称累积误差 向前欧拉公式Vn+=yn+b(xn,yn) 讨论计算一步出现的误差 y=(x)+(x1,(x)=yx)+hy(x) 假定yx),估计m的误差:mm =y(x,)-y=2yx,)+0O(b)=O0) 局部截断误差 局部截断误差主项为”,y(xn) 学学实 (大学数学实验) 误差分析估计欧拉公式的局部断误差 「误算法精度的阶的定义 y(xn)=y(x)+ly(x)+y°(xn)+O(h) 差 个算法的局部截断误差为O(h) 向后欧拉公式yn+=y+hf(xn+t,yn) 析 该算法具有阶精度 局部徽新赣政 r=y(x)-y=-2(x)+O(b)=o(h 向前欧拉公式 O(2) 局部截断误差主项为-5y(x) 向后欧拉公式 O(h2) 梯形向前、向后欧拉公式的平均 梯形公式 y(xn)+O(b)=O(2) 改进欧拉公式O(l) 2阶 经典龙格一库塔公式O(h) 4阶
4 实例1 海上缉私(续) 模型的解析解 2 2 2 2 ( ) ( ) ( ) , ( ) ( ) ( ) c x at y b at y dt dy c x at y b c x dt dx − + − − = − + − − = ( c x ) at y dx dy − − = y at dx dy (c − x) + = 2 2 ( ) d y dt cx a dx dx − = ds b dt = 2 2 ds dx dy = + () () 2 2 2 ( ) 1( ) d y a dy c x dx b dx − =+ dx dy 令 p = b a k p k dx dp (c − x) = 1+ , = 2 p(0) 0 = 实例1 海上缉私(续) 模型的解析解 2 ( ) k 1 p dx dp c − x = + p(0) 0 = 2 1 () c x k p p c − − + += 2 1 () c x k p p c − + −= 1 [( ) ( ) ] 2 cx cx k k p c c − − − = − k ab = / 1 < y(0) 0 = 1 1 2 1 1 [ ( ) ( )] 21 1 1 c c x c x kc k k y kc kc k − − + − = −+ + − − 缉私艇的航线y(x)的解析解 x=c时 缉私艇追上走私船的y坐标 缉私艇追上走私船的时间: 1 2 2 b c t b a = − a=20, b=40, c=15 → t1=0.5 a=35, b=40, c=15 → t1=1.6 2 22 1 kc abc y k ba = = − − 微分方程数值解法的误差分析 数值解法: 计算微分方程精确解y(xn)的近似值yn 按照步长h一步步计算, 每步都有误差; 每一步的误差会逐步积累, 称累积误差. 讨论计算一步出现的误差 假定yn= y(xn) , 估计yn+1的误差: y(xn+1)- yn+1 局部截断误差 误差分析 估计欧拉公式的局部截断误差 y(xn+1)在xn处作Taylor展开: ( ) ( ) 2 ( ) ( ) ( ) 3 2 1 y x O h h y x y x hy x n+ = n + ′ n + ′′ n + 向前欧拉公式 1 (, ) n n nn y y hf x y + = + yn= y(xn) ( ) ( , ( )) ( ) ( ) n 1 n n n n n y = y x + hf x y x = y x + hy′ x + 2 3 2 1 11 ( ) () () () 2 n nn n h T yx y y x Oh Oh + ++ = −= + = ′′ 局部截断误差主项为 ( ) 2 2 n y x h ′′ 误差分析 估计欧拉公式的局部截断误差 ( ) ( ) 2 ( ) ( ) ( ) 3 2 1 y x O h h y x y x hy x n+ = n + ′ n + ′′ n + 向后欧拉公式 ( , ) n+1 = n + n+1 n+1 y y hf x y 2 3 2 1 11 ( ) () () () 2 n nn n h T yx y y x Oh Oh + ++ = − =− + = ′′ 局部截断误差主项为 ( ) 2 2 n y x h − ′′ 梯形 公式 向前、向后欧拉公式的平均 3 4 3 1 () () () 12 n n h T y x Oh Oh + =− + = ′′′ xn xn+1 x y Pn A Q B 误 差 分 析 算法精度的阶的定义 一个算法的局部截断误差为O(hp+1) 该算法具有p阶精度 局部截断误差 精度 向前欧拉公式 O(h2) 1阶 向后欧拉公式 O(h2) 1阶 梯形公式 O(h3) 2阶 改进欧拉公式 O(h3) 2阶 经典龙格-库塔公式 O(h5) 4阶
奥例2翡肉强食 问自然界中同一环境下两个种群之间的生存方式 模型食饵(甲)的密度x(0辅食者(乙的密度 题相互竞争相互依存弱肉强食 x/x=rr>0 甲独立亭的增长率 x/x=r-叫y,a>0己健甲的增长亭魂小 弱肉种群甲靠丰富的自然资源生存食饵(Prey) 视小量与y北万比 强食 y/y=-d,d>0 立存的死亡率d 种群乙靠捕食种群甲为生捕食者( Predator) j/y=-(d-bx),b>0甲宽乙韵死亡小 个种聊的款量如变? 小量与x下比 =(r-ay)x=rx-axy Volterra模型 x(0)=x,y(0)=y x(,y()无解析解 实2翡肉强食型的教解 2弱肉强食型的解析 a=0.1,b=0.02 MATLAB 6.5.1.Ink I-a x(O)=xo,y(0)=yo xo=25,yo=shier.m,shier.m j=-(d-bx)y dy -(d-bx)y 由↓始条确宠 可以相机能是射朗曲能Hx(),y(是屠期函数 证明(c在一定此圆內) 求x()y(0-周期的平均值:x,y j(0=-d-bx)y 4 x(0=(/y+d)/b 猜 测x(0,y(0是期函款;y(x)是封朗曲能 1 Iny(T)-In y(o), d x 数值积分计算一个周期的平均值:x≈25,y≈10 x(odu b (学静学实鉴 (大学数学实验) 实向2翡肉强食桃型的析解 值算法的收煞性和宠性 x(y()一周期的平均值: 收敛性步长h→0时数值解无限接近解析解x 欧拉和龙格库塔方法的共同点只用vn计算ym单步法 r=1d=05,a=01.b=002F=25下=10与计 r~食饵增长率 显式)单步法:{x=+加(x,h 释 结果解 a~捕食者对食饵的捕获能力 单步法有P阶精度局部截断误差O(p+1) x=(r-a ly=-(d-bx)y d-撸食者死亡率 ox,y,h)-(x,,≤少y-(L>0) d单步法收做 b~食饵对捕食者的喂养能力 d整体误差c=xn=Oh)0(b-0) r个a→y个d↑b↓→x个既相工制约 向前欧拉公式改进欧拉公式4阶龙格库塔公式 又相互伉存 燕你溪差en=O(en=O(2) en=O(h)
5 实 例 2 弱肉强食 问 题 自然界中同一环境下两个种群之间的生存方式 相互竞争 相互依存 弱肉强食 弱肉 强食 种群甲靠丰富的自然资源生存 食饵(Prey) 种群乙靠捕食种群甲为生 捕食者(Predator) 两个种群的数量如何演变? 实 例 2 弱肉强食 模型 食饵(甲) 的密度x(t), 捕食者(乙)的密度y(t) x&/ x = r, r > 0 甲独立生存的增长率r x& / x = r − ay, a > 0 乙使甲的增长率减小, 减小量与 y 成正比 y& / y = −d, d > 0 乙独立生存的死亡率d y& / y = −(d − bx), b > 0 甲使乙的死亡率减小, 减小量与 x成正比 0 0 ( ) ( ) (0) , (0) x r ay x rx axy y d bx y dy bxy x xy y ⎧ =− =− ⎪ ⎨ =− − =− + ⎪ ⎩ = = & & Volterra模型 x(t), y(t)无解析解 实例2 弱肉强食 模型的数值解 ⎪ ⎩ ⎪ ⎨ ⎧ = = = − + = − 0 0 x(0) x , y(0) y y dy bxy x rx axy & & 0 0 1, 0.5 0.1, 0.02 25, 2 r d a b x y = = = = = = 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)是封闭曲线 数值积分计算一个周期的平均值： x y ≈ ≈ 25, 10 MATLAB 6.5.1.lnk shier.m, shier1.m d bx y r ay x dy dx ( ) ( ) − − − = 实例2 弱肉强食 模型的解析解 ⎩ ⎨ ⎧ = − − = − y d bx y x r ay x ( ) ( ) & & x e y e c d xx r ay = − − 相轨线 ( )( ) dy y r ay dx x d bx − = − + c由初始条件确定 相轨线是封闭曲线 (c在一定范围内) 求x(t),y(t)一周期的平均值: x , y 可以 证明 y&(t) = −(d −bx)y x(t) = ( y& / y + d ) / b ∫ = T x t dt T x 0 ( ) 1 b d x = x(t), y(t)是周期函数 (周期记作T) ) ln ( ) ln (0) ( 1 b dT b y T y T + − = 实例2 弱肉强食 模型的解析解 r y a x(t),y(t)一周期的平均值: = b d x = r =1, d = 0.5, a = 0.1, b = 0.02 x = 25, y =10 r ~食饵增长率 a ~捕食者对食饵的捕获能力 d ~捕食者死亡率 b ~食饵对捕食者的喂养能力 结果解释 ⎩ ⎨ ⎧ = − − = − y d bx y x r ay x ( ) ( ) & & 与计算结果同 r ↑, a↓⇒y↑ d ↑, b↓⇒x ↑ 既相互制约 又相互依存 欧拉和龙格-库塔方法的共同点:只用yn计算yn+1 1 0 0 ( , ,) ( ) n n nn y y h xyh y yx ⎧ + = + ϕ ⎨ ⎩ = 单步法 步长h→0时数值解yn无限接近解析解y(xn 收敛性 ) 单步法有p阶精度(局部截断误差O(hp+1)) ϕ ϕ ( , , ) ( , , ) ( 0) xyh xyh Ly y L − ≤− > 向前欧拉公式 改进欧拉公式 4阶龙格-库塔公式 数值算法的收敛性和稳定性 (显式)单步法: 整体误差 en=O(h) en=O(h2) en=O(h4) 单步法收敛 整体误差 en=y(xn)-yn=O(hp)→0 (h→0)
大学酸学实 值犷油的收敛位和视定 刚位羽永与刚性方覆 稳定性计算中舍入误差不会随步数的增加无限增大 现象动能或电路系能的款单机盟 的误差&p灬kk=12-日算法穗定 y2=f(xy)中y=(y)+x)+一 dx()=c-m-e-"2+1瞬解与乳 (特征根λ x>0→微分方程稳定 e00快瞬态解计算到产0.005时已衰减到45105 en~慢瞬态解计算到=20时才衰减到45×10 向前欧拉公式“,+机,,=A中=(-边两求稳糟度达到0需算到二20由慢瞬态解2决定 …4.中p-h11dh≤21 态解选取步长h由快瞬态解元=2000定 向后欧拉公式m=男一M日言h任意 龙格-库塔公式h10 求靠则方覆的命令:ode23s,odel5s普(用帽) R(小慢瞬态解mR(A 刚性比 刚性方程 -(10°+1)x1-(10°+2)x2 特征根λmL2106 快瞬态解→→步长充分小 刚性比s=106 慢瞬态解→积分区间长 日传统的数值方法无能为力 x1(0)=10°/4,x2(O)=10°/4-1/2 刚性非线性方程组□心线性常系数方程组 解析解 线性化方法 502+x-→Smm (学静学实鉴 布置实验 用 MATLAB软件掌握求微分方程数值解的方 法,并对结果作初步分析 2通过实例学习用微分方程模型解决简化的实际问 内享课上布量,或参见网络学堂
6 数值算法的收敛性和稳定性 稳定性 计算中舍入误差不会随步数的增加无限增大 yn的误差εn ε n+ k ≤ ε n , k = 1,2,L y′ = f (x, y) λ>0→微分方程稳定 ( , ) n 1 n n n y = y + hf x y 向前欧拉公式 + 向后欧拉公式 nn n 1 1 y y hy + + = − λ 经典龙格-库塔公式 h ≤ 2.785/λ 算法稳定 ( , ) ( , )( ) ( , )( ) * * * * * * * * y f x y f x y x x f x y y y ′= + x − + y − y′ = −λy, λ > 0 (特征根 -λ) n n ε (1 λh)ε +1 = − nk n ε ε + ≤ 1 1 − hλ ≤ h ≤ 2 /λ 1 1 1 n n h ε ε λ + = + h任意 x y ce−λ = n =(1−hλ)y 刚性现象与刚性方程 & x&+kx&+rx= f(t) (k,r >0), x(0) =a, x&(0) =b 现象 振动系统或电路系统的数学模型 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决定 龙格-库塔公式 h10 刚性非线性方程组 线性常系数方程组 传统的数值方法无能为力 线性化方法 刚性方程 刚性现象与刚性方程 刚性方程的MATLAB求解 ode23, ode45 解刚性方程的困难 步长自动变小 计算时间很长 求解刚性方程的命令: ode23s, ode15s 等 (用法相同) 11 2 6 6 2 12 6 6 1 2 2 (10 1) (10 2) (0) 10 / 4, (0) 10 / 4 1/ 2 xx x x xx x x ⎧ = + ⎪ ⎨ =− + − + ⎪ ⎩ = =− & & 例 特征根λ1=-1,λ2=-106 刚性比 s=106 6 6 6 10 1 6 6 10 2 10 ( ) ( 1) 4 10 (10 1) ( ) ( 1) 4 2 t t t t xt e e xt e e − − − − ⎧ =+− ⎪ ⎪ ⎨ ⎪ + =− + + ⎪⎩ 解析解 MATLAB 6.5.1.lnk Stiff.m,stiff1.m 布置实验 目的 1.用MATLAB软件掌握求微分方程数值解的方 法，并对结果作初步分析； 2.通过实例学习用微分方程模型解决简化的实际问 题。 内容 课上布置，或参见网络学堂