《常微分方程课程设计》指导书 5-6

第5章一阶方程组与高阶方程数值方法5.1 一阶常微分方程组数值解方法考虑一阶常微分方程组的初值问题(5.1.1) 若用向量形式表示，可记为，，初始条件，于是(5.1)可写成(5.1.2)(5.1.2)形式上同初值问题(1.2.1)类似，只要看成向量方程即可.因此前面关于单个方程的初值问题数值方法均适用于方程组(5.1.2)，相应理论也可类似地得到.下面仅对如下两个数值方法作说明.梯形法(5.1.3)四阶R-K方法(5.1.4)其中这些公式形式上与单个方程初值问题的数值方法完全一样，但注意这里y及f均为N维向量，在计算机上求解时可直接从数学软件库中选择所要方法，它们都是按方程组初值问题编写的，当N=1时就是前面讨论的单个方程情形.实例演练：自然界中不同种群之间存在着这样一种非常有趣的相互依存、相互制约的生存方式，种群甲靠丰富的自然资源生长，而种群乙靠捕食甲为生，兔子和山猫、落叶松和蚜虫都是这种生存方式的典型.生态学上称种群甲为食饵（Prey ），称种群乙为捕食者（Predator ）,二者共处组成捕食者-食饵生态系统.下面是由意大利数学家V olterra 提出的一个简单的生态学模型.例5.1 食饵甲和捕食者乙在时刻t 的数量分别记作()x t 、()y t ，当甲独立生存时它的（相对）增长率为r ，即()xt rx = ，而乙的存在使甲的增长率减小，设减小的程度与种群数量成正比，于是()x t 满足方程()()xt x r ay rx axy =-=- (5.1.5) 比例系数a 是反映捕食者掠取食饵的能力.捕食者离开食饵无法生存，设乙独自存在时死亡率为d ，即()yt dy =- ，甲为乙提供食物相当于使乙的死亡率降低，并促使其增长.设这个作用与甲的数量成正比，于是()y t 满足方程()()yt y d bx dy bxy =-+=-+ (5.1.6) 比例系数b 是反映食饵对捕食者的供养能力. 设食饵和捕食者的初始数量分别为00(0),(0)x x y y == (5.1.7)微分方程（5.1.5）、（5.1..6）及初始条件（5.1.7）确定了食饵和捕食者数量(),()x t y t 的演变过程. 改写方程组为向量形式，即1112xrx ax x =- 2212xdx bx x =-+ 110220(0),(0)x x x x ==21,[,]xAx A diag r ax d bx ==--+ 1020(0)[,]T x x x =设各参数分别为10201,0.5,0.1,0.02,25,2r d a b x x ======直接用matlab 中的ode45命令求解：function y=shier(t,x) r=1;d=0.5;a=0.1;b=0.02;y=diag([r-a*x(2),-d+b*x(1)])*x; ts=0:0.1:15;x0=[25,2];[t,x]=ode45('shier',ts,x0);plot(t,x),grid,gtext('x1(t)'),gtext('x2(t)'), pause,plot(x(:,1),x(:,2)),grid,xlabel('x1'),ylabel('x2')图5.1 ode45命令求解结果和相轨线例 5.2 '0.01'0.02x x xyy y xy=-⎧⎨=-+⎩0≤ t ≤ 15，x(0)=20, y(0)=20，请用四阶Runge-Kutta 方法来求解，并画出相轨线图.附程序和图形结果显示：f1=inline('x-0.01*x*y'); f2=inline('-y+0.02*x*y'); h=.15;t=0:h:15; x(1)=20;y(1)=20;for n=1:length(t)-1 %方程组的四阶R-K 方法，代公式算之 tn=t(n);xn=x(n);yn=y(n);k11=f1(xn,yn);k21=f2(xn,yn);k12=f1(xn+.5*h*k11,yn+.5*h*k21); k22=f2(xn+.5*h*k11,yn+.5*h*k21); k13=f1(xn+.5*h*k12,yn+.5*h*k22); k23=f2(xn+.5*h*k12,yn+.5*h*k22); k14=f1(xn+h*k13,yn+h*k23); k24=f2(xn+h*k13,yn+h*k23);x(n+1)=xn+h*(k11+2*k12+2*k13+k14)/6; y(n+1)=yn+h*(k21+2*k22+2*k23+k24)/6; endsubplot(121)plot(t,x,'-*b',t,y,'-sr');legend('x(t)食饵','y(t)捕食者'); subplot(122)plot(x,y)；legend('相轨线图')图5.2 四阶R-K法所得数值解曲线和相轨线图5.2 高阶常微分方程数值解方法对于高阶微分方程初值问题，原则上总可归结为一阶方程组，例如下列m阶微分方程(5.2.1)初始条件为(5.2.2)只要引进新变量则可将[WTBX]m阶方程(5.2.1)化为如下一阶方程组(5.2.3)初始条件(5.2.2)则相应化为(5.2.4)实例演练：（导弹追踪问题）设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中？解法一（解析法）假设导弹在t 时刻的位置为P(x(t), y(t))，乙舰位于),1(0t v Q .由于导弹头始终对准乙舰，故此时直线PQ 就是导弹的轨迹曲线弧OP 在点P 处的切线，即有xyt v y --=1'0 即y y x t v +-=')1(0 (5.2.5)又根据题意，弧OP 的长度为AQ 的5倍， 即t v dx y x0025'1=+⎰(5.2.6)由(5.2.5),(5.2.6)消去t 整理得模型:(1)"x y -=初值条件为: 0)0(=y 0)0('=y所求解即为导弹的运行轨迹:245)1(125)1(855654+-+--=x x y 当1=x 时245=y ，即当乙舰航行到点)245,1(处时被导弹击中. 被击中时间为:00245v v y t ==. 若v 0=1, 则在t=0.21处被击中.解法二(数值解)令y 1=y,y 2=y 1’，将方程（3）化为一阶微分方程组122''/(1)y y y x =⎧⎪⎨=-⎪⎩1.建立m-文件eq1.mfunction dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2);dy(2)=1/5*sqrt(1+y(1)^2)/(1-x);2. 取x 0=0，x f =0.9999, [x,y]=ode15s('eq1',[x0 xf],[0 0]);plot(x,y(:,1),’b.')hold ony=0:0.01:2;plot(1,y,’b*')结论: 导弹大致在（1，0.2）处击中乙舰解法三(建立参数方程求数值解)设时刻t 乙舰的坐标为(X(t)，Y(t))，导弹的坐标为(x(t)，y(t)).1．设导弹速度恒为w ，则222)()(w dtdydt dx =+ (5.2.8) 2. 由于弹头始终对准乙舰，故导弹的速度平行于乙舰与导弹头位置的差向量，即⎪⎪⎭⎫ ⎝⎛--=⎪⎪⎪⎪⎭⎫ ⎝⎛y Y x X dt dy dt dx λ， 0>λ (5.2.9) 消去λ得⎪⎪⎩⎪⎪⎨⎧--+-=--+-=)()()()()()(2222y Y y Y x X wdt dy x X y Y x X w dt dx(5.2.10) 3．因乙舰以速度v 0沿直线x=1运动，设v 0=1，则w=5，X=1，Y=t, 因此导弹运动轨迹的参数方程为：⎪⎪⎪⎩⎪⎪⎪⎨⎧==--+-=--+-=0)0(,0)0()()()1(5)1()()1(52222y x y t y t x dtdyx y t x dt dx(5.2.11)4. 解导弹运动轨迹的参数方程 建立m-文件eq2.m 如下： function dy=eq2(t,y) dy=zeros(2,1);dy(1)=5*(1-y(1))/sqrt((1-y(1))^2+(t-y(2))^2); dy(2)=5*(t-y(2))/sqrt((1-y(1))^2+(t-y(2))^2);取t0=0，tf=2，建立主程序chase2.m 如下： [t,y]=ode45('eq2',[0 2],[0 0]); Y=0:0.01:2;plot(1,Y ,'-')， hold on plot(y(:,1),y(:,2),'*')5. 结果见图5.3所示.导弹大致在（1，0.2）处击中乙舰，与前面的结论一致.结论：时刻t=0.21时，导弹在（1，0.21）处击中乙舰。

5.3 练习题1. 对于捕食者-食饵模型'(20 2.50.01)'(100.50.02)u u u vv v v u=+-⎧⎨=--+⎩其中0≤t ≤15，u(0)=20,v(0)=20，请用四阶Runge-Kutta方法和MATLAB命令中的ode45、ode45s来求解，并画出解曲线对比图.2. 对于vdp方程的处置问题2221000(1)0(0)2;'(0)0d x d xx xdt dtx x⎧--+=⎪⎨⎪==⎩请用四阶Runge-Kutta方法求解，并画出解曲线,并和ode45s的解曲线对比.第6章总结1.首先，给出了课程设计的意义和主要内容，同学们学习的重点和实验作业要求.2.本书数值求解问题主要是一阶微分方程及方程组的初值问题(1.2.1)，它是在解存在唯一的前提下讨论的，即对y满足Lips的条件（1.2.2）.求数值解先要建立求解的差分公式就是方程（1.2.1）进行离散化.涉及的基本概念是方法的局部截断误差及其主项，方法的阶，收敛性与整体误差，和方法的绝对稳定性.得到的方法按使用初值个数不同分为单步法及多步法两类.单步法计算时只要用到前一点的值，多步法则要用到前面的多个值.另外按是否用到的值方法又分为显式方法和隐式方法，对显式方法计算较为方便，而对用到的隐式方法计算时还要利用迭代，当然实际计算时多采用预测－校正公式.3.本文重点是根据局部截断误差定义建立单步法及多步法的公式和求出其局部截断误差，并确认公式的阶.主要方法是利用Taylor展开，对二阶以上的显式单步法，即显式R-K方法用二元函数的Taylor公式（3.2.14），其余所有的公式都可包含在多步法（K＝1为单步法）中，都可用一元函数的Taylor展开方法求得，其方法是将公式局部截断误差的右端各项在处按Taylor展开并令，则可得到公式系数满足的方程组，从而可求得所要建立的数值方法公式及其余项.4.本文构造的简单方法有Euler法，后退Euler法，梯形法、单步法有二阶，三阶，四阶的R-K方法以及多步法Adams显式与隐式方法，Milne法与Hamming 法等，要能用这些方法计算给定方程的数值解，对单步法要能给出它们绝对稳定区间，并应用于步长h的选取，对隐式方法要能用收敛的迭代法计算出解，也能用预测-校正方法求解.参考文献[1]周义仓，靳祯，秦军林. 常微分方程及其应用[M]. 北京：科学出版社，2003.7[2]徐萃微，孙绳武. 计算方法引论(第二版)[M]. 北京：高等教育出版社，2002.1[3]封建湖. 计算方法典型题分析解集[M].西安：西北工业大学出版社，2000[4]冯康. 数值计算方法[M]. 杭州：浙江大学出版社，2003[5]李庆扬，关治，白峰杉. 数值计算原理.北京：清华大学出版社，2000[6]任玉杰. 数值分析及其MATLAB实现[M]. 北京：高等教育出版社，2007[7]李庆扬，王能超，易大义. 数值分析[M].北京：清华大学出版社，施普林格出版社，2001[8]胡健伟，汤怀民. 微分方程数值解法（第二版）[M]. 北京：科学出版社，2007[9]李立康，於崇华，朱政华. 微分方程数值解法[M]. 上海：复旦大学出版社，1999[10]C.W.吉尔. 常微分方程初值问题的数值解法[M]. 费景高译.北京：科学出版社，1978[11]马知恩，周义仓. 常微分方程定性与稳定性方法.北京：科学出版社，2007.1附注：网易邮箱mathofhdu_angels@为公共信箱密码为：mathofhdu。