第一章等切面曲线和相似曲线 W. gander. s. Barton and. hrebicek 1.1引言 在这一章里,我们将用 MATLAB解两个相似的微分方程系统.首先,我们推广古典的等切面 曲线问题,计算小孩拉玩具的轨迹.然后计算狗攻击慢跑者的轨迹.我们还用 MATLAB显示这些 运动轨迹 12古典等切面曲线 在17世纪, Gottfried wilhelm leibniz曾讨论下面问题,参见12,1:将一个手表系在一条链 子上,当沿直线拉此链的一端时,在平面上,此手表所描述的轨迹是什么? 假设a为此链的长度,手表看成一个点,如果假设手表开始位于r轴的点(a,0)处,我们从 原点开始在y轴的正方向上拉,则容易求解此问题(2],(参见图11) 图1.1古典等切面曲线 手表 从图1.1,我们可得函数y(x)的微分方程 √a2-x2 为了解方程(11),我们只需积分: y:=-int(sqrt(a"2-x"2)/x, x)+c
W. Gander, s. Barton and J. hfebicek cr2+a arctan( (1.2) 当执行不定积分时, MAPLE不包含常数.因此我们增加一个积分常数c.利用 lindsay(x)=0.我 们能确定它的值 C:= solve(limit(y, x=a),c); e:=Ⅰa-丌 利用 MAPLE53,我们从int直接得到一个实表达式,因为对于实表达式,此常数为0,我们不必 计算一个复的积分常数 yold: =-sqrt(a"2-x2)+a*arctan(sqrt(a"2-x"2/a) yold: =-va- 2-22+a'arctanh( (1.3 为了证明两个结论相同,我们把它们相减并转成指数对数形式,化简并展开 e:=expand(simplify(convert((y-yold)/a, expln))); e:=-ln(√a2-x2-a)+1x+1m(√a2-x2+a2) 在 MAPLE中,不可能再化简此式,但是应该注意,第一个对数函数中的量是第二个的负值.因此 如果用一个正的未知数d代替a-√a2-x2,类似地用一d代替va2-x2-a,则 MAPLE能化简这 个表达式: > simplify(subs({(a2x-2)^(1/2)-a=-d,-(a^2-x-2)(1/2)+a=d},e)); 假设被拉物体的初始位置在y轴上的点(O,a),我们从原点开始拉,但这次在x轴的正方向上 考虑在物体的轨迹上的点(x,y(x)在x轴上,此链的端点是(x-y(x)/y/(x),0),即切线与x 轴的交点(与 Newton迭代一步所得的点相同!.所以,具有常数长度a的此链,导出微分方程 v(x)2 不能通过求积分直接解它.因此,需要调用微分方程的求解程序 dsolve, >assume(a>0) eq (y(x)/diff(y( (a y(a y(x)2
第一章等切面曲线和相似曲线 3 p: dsolve(eq,y(x)): p:=V-y(x)2+a-2-aarctanh( +x=C1 y(x)2+ y(e)2+a"+a arctan( y(x)2+ 并得到两个解.从初始条件y(0)=a和物理学,对于a>0,可得y(m)>0和y(x)-_C1: limit(lhs(subs(x=o, y(x=y), P[1])),y=a) 于是,解满足方程 y(a)"+a-a arctan( )+x=Ia^丌 V-y(x)2+ 由 MAPLE53,可得到实表达式 a2-y(a)"-a arctan 当然,我们也可在方程(12)和(1.3)中,互换变量和y得到这些方程.值得注意的是,因为对 x=0,有奇异性y(0)=∞,所以用数值方法解方程(14)是困难的 13小孩和玩具 让我们解决一个更一般的问题,假设一个小孩在平面上沿一曲线行走,此曲线由两个时间的 函数X(t)和Y(t)确定 假设此小孩借助长度为a的硬棒,拉或推某玩具,当此小孩沿曲线行走时,我们计算玩具的 轨迹.设(x(t),y(t)是玩具的位置 从图1.2可得下列方程 1.(X(t),Y(t)与(x(t),y(t)之间的距离总是硬棒的长度,于是 (X-x)2+(Y-y)2 (1.5) 2.玩具总是在硬棒的方向上运动,因此,两个位置的差向量是玩具的速度向量的倍数,vr X-a 其中入>0. (16)
W. gander s. Barton and. Hrebicek 图12速度vc和vT [x(),y(I)] 小孩的速度v 玩具的速度Ⅵ 3.玩具的速度依赖于小孩的速度向量vc的方向,例如,假设小孩在半径为a(硬棒的长)的圆 上行走.在此特殊情况下,玩具停留在此圆的圆心,根本不运动(这是第一个例子的最后状 态,参见图1.3). 从图1.2可知,小孩的速度vω在硬棒上的投影的模是玩具的速度ⅵr的模 将方程(1.6)代入方程(15)中,可得 √a2+y2 于是 (1.7) 为了得到x和,我们要解方程(1.7).因为玩具的速度的模|vr|=|v。 cosa,见图1.2,这由下面步 骤可得到 ·标准化差向量(X-x,Y-y),可得单位长的向量w 确定vC=(x,Y)在w生成的子空间上的投影.因为vw=| vcllw cos a和wl=1,所以 这就是内积vw 现在,我们能用 MATLAB写出函数,求解微分方程系统 算法1.1函数f function zs =f(t, z) LX Xs Y Ys]=child(t) ⅴ=[Ks;Ys] w=[X-z(1);Y-z(2)] 函数f调用一个函数 child,对于给定的时间t,它返回小孩的位置(X(t),Y(t)和速度(Xs(t),Ys(t) 例如,考虑小孩在圆X()=5cost;Y(t)=5sint上行走,此时相应的函数chi1d是 算法1.2函数 Child function [X, Xs, Y, Ys]= child(t)
第一章等切面曲线和相似曲线 5 5*sin(t) Xs =-5*sin(t): Ys 5*cos(t) MATLAB提供两个M文件ade23和ode45,用于求积微分方程.在下面的主程序中,我们将调用 这些函数之一,并定义初始条件(注意,当t=0时,小孩在点(5,0)处和玩具在点(10,0)处) >>% main1 >>y0=[10o] >>[ty]=ode45(f’,[0100,yo) >> clf hold on; >>axis([-610-610]); >>axis >>plot(y(:,1),y(:,2)) 图1.3小孩在圆上行走 如果绘制两列y,则可得玩具的轨迹(参见图13),而且在同一图中,用语句加进小孩的曲线: >>t=0:0.05:6.3; >>[X, Xs, Y, Ys]= child(t) >>hold off 注意,在程序中,硬棒的长a没有明显出现;它由玩具的位置(初始条件),和t=0时小孩的位置(函 数chi1d)隐含地定义 我们用几个例子来结束本节.假设小孩沿sin函数的图形行走:X(t)=t和Y(t)=5sint,小 孩的曲线用虚线表示,初始条件c(0)=0和y(0)=10,我们得到图14 另一个例子,小孩在圆X(t)=5cost,Y(t)=5sint上行走,初始条件x(0)=0和y()=10, 我们得到一个象花的玩具轨迹(参见图1.5)
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 > 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)))
