W. Gander. S. Barton and J. Hrebicek 图1.4例2 1.4慢跑者和狗 我们考虑下面问题:为了每天锻炼,一个慢跑者在平面上沿着他喜欢的路径跑步,突然一只 狗攻击他,这只狗以恒定速率w跑向慢跑者,计算狗的轨迹 狗的轨迹具有如下性质:在任意时刻,狗的速度向量都指向它的目标慢跑者.假设慢跑者在 某路径上跑步,他的运动由两个函数X(t)和Y(t)描述 假设当t=0时,狗是在点(x0,o)处,在时刻t时,它的位置是(x(t),y(t),下列方程成立 2+y2=v2:狗以恒定速率跑 2.狗的速度向量平行于慢跑者与狗的位置的差向量 X 其中λ>0. 如果将它代入第一个方程,我们得到 (x=) 求解此方程,可得A: 最后,将此式代入第二个方程,可得狗的轨迹的微分方程:
第一章等切面曲线和相似曲线 花形轨迹 15-10-50 我们将利用M文件ode23.m或ode45.m,求积微分方程系统.注意,当狗达到慢跑者时,系统(1.8) 有一个奇异点.此时差向量的范数变为0,我们必须停止积分.上面提到的求积微分方程的 MATLAB 函数,要求输入独立变量的一个区间.现在, MATLAB5.0也可定义另一个求积终止准则,它不同 于给定独立变量的一个上界,而是通过检查函数的零交点,来终止积分.在我们的例子中,当狗达 到慢跑者,即‖(X-x,Y-y川‖变成很小时,将终止积分.为此,在M函数dogm中,我们必须 增加第三个输入参数和两个新的输出参数.积分器ode23或ode45,以两种方法调用此函数:第 种方法是取消第三个参数,函数只返回参数zs:狗的速率.第二种方法是,将关键字‘ events赋 值给参数红1ag,此关键字告诉函数,在第一个输出zs中返回零交点函数.第二个输出 asterina1 是一个逻辑向量,它告诉积分器,当第一个输出的分量变成0时,它们迫使过程终止.具有这种 性质的每一个分量,在 asterina1中,用非零标记.第三个输出参数 direction也是一个向量, 它表示:对z8的每个分量,零交点是否仅被视为增加值( direction=1),下降值( direction 1)或两者( direction=0)零交点的条件将在积分器里检查,dog和主程序中,必须声明狗的 速率W是全局变量.M函数 Jogger,m给出慢跑者的轨迹 算法13函数Dog function [zs, isterminal, direction] =dog(t, z, flag); global w %w= speed of the dog nh= norm(h);
8 W. Gander. S. Barton and. Hrebicek if nargin <31 isempty(flag)% normal output zs=(w/nh)*h 18e switch(flag) case 'events, at norm(h)=0 there is a singularity -1e-3;% A this is a stopping event direction=o: don't care if decrease or increase otherwise error (['Unknown flag: flag]); end 主程序main2.m定义初始条件,调用ode23进行积分.为了求积,我们必须提供时间t的上界 >> main2.m >> global w >>yo=[60: 70]: %initial condition, starting point of the dog >>W=10; %w speed of the dog >> options= odeset('RelTol,1e-5,' ' Events','on') >>[t, Y]= ode23('dog, [o, 20], yo, options) > clf: hold o >>axis([-10,100,-10,701); >>plot(Y(:1),Y(:,2)); >>J=[] >>for h=1: length(t), jogger(t(h)) J=[J;w’]; >plot(J(:,1),J(:2),:); 如果达到时间t的上界或狗赶上慢跑者,积分将终止.对于后者,我们设置ODE选项的标志 Events 为‘on’,这告诉积分器检查,具有红1ag= events,的被调用函数dog的零交点.调用ode23后, 变量y包含有两个函数x(t)和y(t)的值的表.用语句plot(Y(:,1),Y(:,2)),我们画出狗的轨 迹.为了显示慢跑者的轨迹,我们需要用向量t和函数 Jogger计算它 现在,我们计算几个例子.首先,假设慢跑者沿x轴跑步 算法1.4第一个慢跑者例子 function s= jogger(t) s=[8*t;O] 在上面的主程序里,我们选取狗的速率为=10,因而X(t)=8t,慢跑者是较慢的.如图16所 狗抓住慢跑者
第一章等切面曲线和相似曲线 如果想指出慢跑者有麻烦的位置,(也许设立一个小纪念物,)我们利用下面文件 cross.m 算法1.5画十字标 function cross(Cx, Cy, v) l dra Cx. C ross of height 2.5 and width 2*v Kx = [Cx Cx Cx Cx-v Cx+v] Ky =[CyCy+2. 5*V Cy+1. 5*v Cy+1.5*V Cy+1. 5*v] 图中的十字标是由主程序增加语句 >>p= max(size(y)) s(Y(P,1),Y(p,2),2) > hold off: 所产生.下一个例子显示慢跑者转向的位置并跑回家 算法1.6第二个慢跑者例子 unct1。ns agger ift6,s=[8*t;0] [8*(12-t);0 图16慢跑者沿直线y=0跑
10 W. Gander, S. Bartori and. hrebicek 图17慢跑者向回跑 用前面相同的主程序,当t=93时,狗赶上慢跑者(参见图17)现在,考虑一个更快的慢跑 者在一个椭圆上跑步 算法17第三个慢跑者例子 function s=jogger(t): s=[10+20*cos(t) 如果狗跑得快(u=19),当t=8.97时,它赶上慢跑者(参见图1.8).最后,考虑一只老的、慢的狗 (=10),它在一个椭圆上跑,并且努力赶上慢跑者.然而,狗不在椭园上的某点等待,而是在它 的目标后面跑(太慢),我们能看到一个稳定的状态,狗在椭圆内一个闭轨迹上跑(参见图19) 15用 MATLAB显示运动 同时显示小孩和玩具,或狗和慢跑者的运动,而不只是画出它们静态的轨迹,这样效果会更 好,这可用 MATLAB的图形句柄命令.关于小孩和玩具的主程序如下 > options= odeset('RelTol', 1e-10) >>[t y]= ode45('f,, [o 40], yo, options); >>[X, Xs, Y, Ys]=child (t) xmin =min (min (X), min (y (: 1)))