F. Klvaria init =x(o) dx(o)=vo*cos(theta), dz(o)=vo*sin(theta): equations of motion of ball in air rotation (drag force only) eqnto: diff(x(t),t)=dx(t) diff(dx(t),t)=-0508*alphadx(t)*v(t) iff((t),t)=dz(t) diff(dz(七),七)=-g-0.508*a1pha*dz(t)*v(t) equations of motion of rotating ball in air # (influence of Magnus effect) eqnti: diff(x(t),t)=dx(t), diff(dx(t), t)=(-Cd*dx(t)+ Cm*dz(t))*alpha*v(t) diff(z(t),t)=dz(t) diff(dzt), t)=-g(Cd*dz (t)+ Cm*dx(t))*alpha*v(t) #numeric values of initial parameters h: =1: vO 25: theta Pi/180*15: theta= 15 degrees solution for non rotating ball reso = dsolvedleqnto, init, var, numeric); 0:=proc(rkf45x)…end > result is a set of equations >res0(0.5); t=.5,x(t)=10.76902243898168,z(t)=2745476721738618, dx(t)=1931517542170,dz(t)=.758480487012169] solution for rotating ball res1: =dsolve(eqnt1, initc), var, numeric) proc(rkf45-r) 函数 dsolve(., numerIc)返回带一个参数(独立变量)的函数.调用此函数得到一组方程 (在此情况下,变量是t,x(t),x(t),vx(t),U-(t).为了访问数值,我们仍用函数subs 绘制图形 作为解法的最后部分,我们将在同一图中,画出三个模型的轨道.为此我们必须解一个求飞 行时间(方程x(t)=0的解)的子问题 为了计算真空中球的飞行时间 marid,我们用(26)的符号解resd和函数 fsolve calculation of tmax -time of falling to earth for ball in vacuum tmaxid fsolve(subs resi tmad:=1,458903640
第二章旋转网球的轨迹 在其它情况下,因为通过对微分方程积分,可得到z的导数U2,所以我们容易用 Newton法求解方 程z(t)=0.于是方程x(t)=0的近似解的递推关系是 算法2.1函数zzex zzero := proc (u, to, z, dz) local t, ts, up using diff(z, t)= subs(u(t),dz) hile abs((tn -ts)/tn)>10"(-4)do; tn: =ts -subs(up, z)/subs(up, dz) 我们可用 marid作为初始近似解to. MAPLE函数 zero实现此 Newton迭代(参见算法21) calculation of the flight time for the other model >tmax=zero(reso, tmaxid, z(t),dzt)); tmax1: =zero (res1, tmaxid, z(t),dz(t)) mac0:=1.362013689 ma1:=9472277855 为了从数值解生成空气中球的轨道,最简单(可能最快)的方法是用一个数组[a(t),x(t1)作为函数 P1ot的输入参数.对于一个具有常数时间步长的区间,为了从 dsolve(.…, numeric)的结果得 到此数组,我们定义简单函数 tabler.为了产生红,蓝,黑色的图形,剩余的 MAPLE语旬是 making graphs Gid: =plot([subs (resid, x(t)), subs (resid, z(t)) t=0.. tmaxid], linestyle = 2) for models with ric solution s [x(t), z(t)] for plotting tablepar proc (u xmin, xmax, npoints) local i, Step; [seq([subs(u(xmin i*Step),x),subs(u(xmin+ i*Step),y)] 1=0.. npoints) Go plot(tablepar(reso, x(t), z(t),o, tmax, 15) >G1: plot(tablepar(res1, x(t), z(t),o, tmax1, 15)
plotting of all graphs plots [display](Gid, GO, G1) 用 MATLAB计算轨道,可得到相同的这些图形(参见图22 算法2.2必要的M函数集合 function xdot= tennisip(t, x, flag) [x(3) g isop(t, x) alpha x(4 -alpha*O 508*x(3)*v g-a1pha*0.508*x(4)*v]; function dot= tennis1p(t, x) cd=(0.508+1/(22.503+4.196*(v/w)"0.4))*a1pha*v cm=etha*/(2.022*+0,981*v)*a1pha*; g-Cd*x(4)-Cm*x(3) 2.3 MATLAB解法 由于问题是非线性的,我们必须用数值方法求解.象 MATLAB这样的数值系统比符号系统更 合适.因此我们将用 MATLAB求解问题 为了求解n个微分方程的系统的初值问题 d ∫(t,c) (210) 其中z=(x1,x2,…,xn,我们用 MATLAB函数ode23(或ode45),实现2阶和3阶(或4阶和5阶) 的嵌入式 Runge-Kutta方法
第二章旋转网球的轨迹 我们必须对每个模型定义M函数,以实现微分方程系统(210).为了在M文件之间传递参数 g,a,w和η,在我们的程序中,定义它们为全局变量.用如下映射 x→c(1),z→x(2),vx→x(3),tx→x(4) 所需的M文件在算法(22)中 在下面主程序中,我们计算和绘制所有三个模型的轨道(参见图22).为了画出球在空气中的 轨道,计算一个充分密集的表{x,,我们用z=z(x)的100个解点作样条函数插值,所用函数为 SP A Trajectory of spinning tennis ball %主nit主a1 ization > global g alpha w eta Z basic constants in MKS units >g=9.81;d=0.063;m=0.05;rho=1.29 >>alpha= pi*d"2/(8*m)*rho > eta 1 >>w=20 A initial conditi >>h=1;v0=25; theta=pi/180*15 >>xin =[o, h, vo*cos(theta), vo*sin(theta)] %±1ig >> tmaxid =(xin(4)+ sqrt(xin (4)2+2*g*xin(2)))/g; A solution in vacuum >>[tid, xid] ode 23('tennisip, [o tmaxid], xin); > [to, xo]= ode23('tennisop, [o tmaxid], xin) soluti >>[t1, x1]= ode23('tennis1p, to tmaxid], xin); >I=max(xid(:,1));x=0:N/100:N (to, max(xid(:, 1)),0, max(xid(:, 2))]) >> hold on >>plot(x, spline(xid(:, 1), xid (: 2),x),2:r) >p1ot(x, spline(xo(:,1),x0(:,2),x),,--b) >p1ot(x,sp1ine(x1(:,1),x1(:;2),x),,-w3) >> hold off 注意,在空气中的两个模型,我们不必计算飞行时间.对于最大的轨道(在真空中的球),利用紧跟 axis语句后的hold语句,轨道被截断,并且不画出x轴以下的图形 如果希望计算飞行时间,类似于 MAPLE解法,我们可以求解.因为从开始就必须对微分方程 组数值积分,得到函数值和导数,这种解法要花费很多的计算时间 然而,用逆插法计算飞行时间的近似值更简单.为此我们要找一个时间区间,在此区间上函 数是可逆的.我们取表中的某些值,在这些值处,轨道的z值是单调递减的.这可由函数min和
F. Klara max 4 Determine flight time by inverse interpolation [z,j=min(x0(:;2);[y,i=ax(x0(:,2)) >>tmax0=sp1ine(x0(i+1:j,2),t0(i+1:j),0) >[z,j-in(x1(:,2));[y,i=ax(x1(:,2)) >> tmax1 spline(x1(i+1: j, 2), t1(i+1: j), o) 为了使区域z(i+1:j)是单调的,我们取第(a+1)个点为z的单调区域的第一个点.于是得到飞 行时间的近似值tmax0=1.3620和 marl=0.9193,比 MAPLE的结果好 图22 2.4 MATLAB5更简单的解法 MATLAB5允许对微分方程求积,直到所得的轨道交于x轴2,因此不需要样条插值.根据第 章第7页的解释,为了检查零交点,必须调整算法(22)中定义的M函数;算法(23)给出修改 的M函数.计算三个轨迹的主程序变得更短:选项 Events设置成on,使积分器检查z轴的零交 点,如果有一个,则停止积分.因此可由积分器返回的时间向量t和t1的最后部分,得到轨迹的 时间. 算法23在 MATLAB5下重写上一节的M文件 function [xdot, isterminal, direction]= tennisip2(t, x, flag) f nargin <3 I isempty(flag)% normal output einhard Jaschke给出此解