第一章等切面曲线和相似曲线 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)<0.因此,第 一个解是我们问题的正确答案 P[1] V-y(a)2+a-a arctan( )+x=C1 我们得到隐式的解v(x).为了求积分常数c1,我们利用约束lmx→0(x)=a >-_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)