Matlab软件与基础数学实验程序(部分)

追击问题：

一敌舰在某海域内以椭圆轨迹航行，其在时间t时刻的坐标为：

x(t)=10+20cost

y(t)=20+5sint

我方战舰恰位于原点处，我战舰向敌舰发射制导鱼雷，鱼雷的速率为20，其运行方向始终指向敌舰，试问敌舰航行在何处将被击中？

2. 若敌舰的运行轨迹变为

x(t)=10+20cost

y(t)=20+20sint

试问敌舰航行在何处将被击中？(无法击中)

3. 若敌舰的运行轨迹变为

x(t)=10+20cost

y(t)=20+20sint

鱼雷速率提高至21，结果如何？

%Matlab程序：

clear;clc

h=;%时间步长

k=1;

t(1)=0;x(1)=0;y(1)=0;%初始值

r=10;

while r>= % k<=250 %

m=(20+5*sin(t(k))-y(k))/(10+20*cos(t(k))-x(k) ++;

if 10+20*cos(t(k))-x(k)>=0

x(k+1)=x(k)+20*h/sqrt(1+m^2);

else

x(k+1)=x(k)-20*h/sqrt(1+m^2);

end

if 20+5*sin(t(k))-y(k)>=0

y(k+1)=y(k)+20*h/sqrt(1+1/m/m);

else

y(k+1)=y(k)-20*h/sqrt(1+1/m/m);

end

r=(x(k)-10-20*cos(t(k)))^2+(y(k)-20-5*sin(t(k)))^2; r=sqrt(r);

t(k+1)=h*k;

k=k+1;

plot(10+20*cos(t(k)),20+5*sin(t(k)), 'r*')

hold on

axis([-10 32 -3 30]);

plot(x,y, 'o')

pause

end

t=t(end),x=x(end),y=y(end)

t =

x =

y =

第二问：速度相同无法击中

第三问：

t = x = y =

%Matlab程序：

clear;clc

h=;%时间步长

k=1;

t(1)=0;x(1)=0;y(1)=0;%初始值

r=10;

while r>= % k<=250 %

m=(20+20*sin(t(k))-y(k))/(10+20*cos(t(k))-x(k) ++;

if 10+20*cos(t(k))-x(k)>=0

x(k+1)=x(k)+22*h/sqrt(1+m^2);

else

x(k+1)=x(k)-22*h/sqrt(1+m^2);

end

if 20+20*sin(t(k))-y(k)>=0

y(k+1)=y(k)+22*h/sqrt(1+1/m/m);

else

y(k+1)=y(k)-22*h/sqrt(1+1/m/m);

end

r=(x(k)-10-20*cos(t(k)))^2+(y(k)-20-20*sin(t(k)))^2; r=sqrt(r);

t(k+1)=h*k;

k=k+1;

plot(10+20*cos(t(k)),20+20*sin(t(k)), 'r*')

hold on

axis([-12 32 -2 42]);

plot(x,y, 'o')

pause

end

t=t(end),x=x(end),y=y(end)

课本P81

1. 某农夫有一个半径10米的圆形牛栏，长满了草.他要将一头牛栓在牛栏边界的栏桩上，但只让牛吃到一半草，问栓牛鼻的绳子应为多长？

设拴牛的绳子长为r, 以圆形牛栏C1 的圆心为原点建立直角坐标系, 见图1, 不妨设拴牛的栏桩为图1中圆形牛栏C1 上的B 点, 其坐标为(10,0), 则所求问题转化为: 求出r,使得以B 点为圆心, 半径为r 的圆C2 与圆C1 相交部分的面积 是圆C1 面积的一半。

解法一：

由于圆形牛栏C1 和圆C2 的方程分别为:

C1: x 2+y 2=100

C2: (x-10)2+y 2=r 2 (1)

联立方程C1, C2, 可得两交点分别为：

⎪⎭⎫ ⎝⎛--2240020,20110r r r C ，⎪⎭

⎫ ⎝⎛---2240020,20110r r r D 设牛吃草的面积为S ，即圆C1与C2的相交部分，则根据题意，S 应为圆C1面积的一半，即

()

2102

1⨯=πS 由图可知，S 的面积可由下面的定积分计算得到： ()212I I S +=

⎪⎪⎭⎫ ⎝⎛-+--=⎰⎰---dx x dx x r r r r 102010222010102222

10)10(2 对于上式中的积分1I ，令u=x-10，则上式可化简为 S ⎪⎪⎭⎫ ⎝⎛-+-=⎰⎰---dx x dx u r r r r 10201022202222

102 上式通过简单积分运算可化为

f=inline('-r^3/800*sqrt(400-r^2)*r^2*asin(r/20)+*pi*r^*sqrt(400*r^2-r^4)+r^2/800*s qrt(400*r^2-r^4)-50*asin(1-r^2/200)','r');

fzero(f,10)

ans =

另解：

S0=1/2*pi*10^2;

r=10; %给定r 的迭代初值;

S=314;

while abs(S-S0)>

r=r+; %r 的值增加;

u=- r::- r^2/20;

s1=trapz(u,sqrt(r^2- u.^2));

x=10- r^2/20::10;

s2=trapz(x,sqrt(10^2- x.^2));

S=2*(s1+s2);

end %当误差小于eps 时, 循环结束;

error=abs(S0-S) %显示误差;

[r,S] %显示r 的值和面积;

error =