第二章旋转网球的轨迹 xdot=[x(3),x(4),0,-g]’; case 'events' at I Ill=o there is a singularity otherwi flag function [xdot, isterminal, direction]= tennisOp2(t, x, flag) if nargin <3 I isempty(flag)% normal output v=sqrt(x(3)2+x(4)^2); xdot=[x(3),x(4), a1pha*0.508*x(3)*V,-g- alpha*o.508*x(4)*v]; switch(fla case 'events,% at I lhll=0 there is a singularity isterminal= 1 eror([ Unknowm f1ag’’f1ag end function [xdot, isterminal, direction]= tennis1p2(t, x, flag) global g alpha if nargin 3 I isempty(flag)% normal output v=sqrt(x(3)2+x(4)^2); cd=(o508+1/(22.503+4.196*(v/v)^0.4))*a1pha*v c*x(3)+cm*x(4),-g-cd*x(4)-Cm*x(3)]); switch(flag) case 'events,% at I l=o there is a singularity
F. Kwania = direction= error(['Unknown flag ,),flag 777. ']); end 积分器取相当大的积分步长,使得轨迹的图形没有样条插值所得的光滑.为了使这些轨迹光滑, 我们必须取小步长.用选项 Maxstep规定步长的上界来实现.通常设此值为时间区间长的十分之 在下列算法中,我们令它为 tmaxid/100,可得到如图22所示的同样光滑的轨迹 f Trajectory of rotating tennis ball h initialization > global g alpha w etha g=9.81;d=0.063; 0.05;rho=1.29 > alpha pi*d"2/(8*m)*rho; >> etha =1 Z initial conditions >>h=1;vo=25; theta=pi/180*15 >>xin [o, h, vO*cos(theta), vo*sin(theta) h flight time for vacuum >>tmaxid =(xin(4)+sqrt(xin(4)2+2*g*xin(2)))/g A setting options for the integrator >>optiona=odeset('Events,'on','Maxstep,tmaxid/100) >> >>[tid, xid]=ode23('tennisip2', [o tmaxid],xin, options) h soluti >>[to, xo] = ode23('tennisop2', [o tmaxid], xin, options) A solution with spin >>[t1, x1]= ode23('tennis1p2', [o tmaxid], xin, options) >c1f; >>axis([0,max(xid(:,1))*1.13,0,max(xid(:,2))*1.13]); >>hold on: > plot(xid(:, 1), xid(:, 2) >>plot(x0(:,1),x0 >>plot(x1(:,1),x1( -g') >>hold off Z flight time for the trajectory without spin >>tmax to(length(to)) h flight time for the trajectory with spin
第二章旋转网球的轨迹 >>tmax1 =t1(length(t1)) 参考文献 [1] E.G. RICHARDSON, Dynamics of Real Fluids, Edward Arnold,1961 [2]A. STEPANEK, The Aerodynamics of Tennis Balls- The Topspin Lob, American Journal of Physics, 1988,pp.138-14
第三章道路照明问题 S Barton, D Gruntz 31引言 这一章我们研究由两盏路灯照明的一条水平的道路,其中P是路灯的亮度,hz是灯的高度 两盏路灯的坐标分别是(0,h1)和(s,h2),其中s是两灯之间的水平距离.令X=(x,0)是两灯之 间道路上的一个点.这一章我们将找出具有最小照明强度的点X.图31给出了在这一章将要讨 论的问题的一个示意图 图31道路照明问题示意图 在第一节我们将在给定两灯的高度和亮度的前提下求出X点的位置.第二节我们通过改变第 二盏灯的高度以极大化X点的照明强度.在最后一节我们还将进一步优化与两灯高度有关的X 点的照明问题.事实上,我们将通过改变两灯的高度来优化道路照明问题 32道路上的最低照明强度的点 这一节我们将寻找两个光源之间照明强度最低的点X.由物理学的知识可知:被光线照射的 物体的亮度依赖于它与光源之间的距离平方的倒数和光线的投射角度,参见[1,2 从X到第一个光源的水平距离是x,到第二个的水平距离是s-x.利用勾股定理我们可以 确定点X到两个光源之间的距离r1和 r=h2+x2,r2=h2+(s-x)2
第三章道路照明问题 从两盏灯到X点的光线强度分别为 P2 1)=n=好+2,2()=n=+(5-x 如果光线的投射角度分别为α1和α2,道路的照明情况还依赖于sina和sina2 h2 于是X点的总的照明强度C()是 C(=)=I1(e)sin a1 + I2(r)sin a= Pahs (31) vG+3√(+(-x) 对于0≤x≤s,极小化C(x)就得到X点的坐标.为了求出极小值,可以对C(x)求导并求出它的 根.我们尝试使用 MAPLE来做此事 >s[1]:=P[1]*h[1]/(h[1^2+x2)-(3/2) >S[2]:=P[2]*h[2]/(h[2]2+(s-x)-2)-(3/2) >C:=S[1]+S[2] =diff(c, x) dc. a P1 hi r 3P2h2(-28+2x) )5/22(h2+(-x)2)5 solve(dc=0,x) 如果你试图执行这个指令,你将会看到MAPE决不会给你回答.使用代数学的处理方法,方程 dC=0能够被转换为关于x的多项式方程.我们把方程dC=0的一项移到方程的右端,两边平 方,再把右端移回左端,通分.则分子必须等于零 diff(s[1], x)"2- diff(s[2], x)"2 eq:=9 (h12+x2)54(h2+(-x)2) collect (prim 最后一个指令的结果是一个关于x的12次的多项式 (P2h2-P2h2)x2+(2P2h2s-10P1hs)x1+…-P2h2h1°。2=0. (32) 式中的常数不确定时,这个多项式很难解出其封闭的形式 我们考虑如下的数值:P1=2000,P2=3000],h1=5m],h2=6m],s=20m]函数 C(x),C(x)≡dC,S1(x)和S2(x)(每一盏灯分别投射到道路上的照明强度)可以使用 MAPLE在图 上绘制出来(见图3.2).从图上可以辨别出包含函数C(x)的零点的区间.我们将在这个区间上用 fsolve求出X的数值. P[1]:=2000:P[2]:=3000:s:=20:h[1]:=5:h[2]:=6 plot({c,S[1,S[2],d},x=-8/2..8*3/2) Xm : fsolve(dc=0, x, 5..10): Xm:=9.338299136