ep3_p7.m

此代码展示了Euler计算方法与差分方法

%  Ep3_p7: << 实验三 >>  敌舰作蛇形逃逸（自补充内容）--- 图形演示版
%  Designed by FGH

h= 0.0014;
H= 120; Vw= 450; Ve= 135;

clf
axis([-5 60 -10 150])
hold on
title('敌舰作蛇形逃逸')
plot(0,H,'bo')
plot(0,0,'r.')
pause

clear Xw Yw Xe Ye;

tk= 0; s= 0;
k= 1;
Xw(1)= 0; Yw(1)= 0;
Xe(1)= 0; Ye(1)= H;
while sqrt( (Xw(k)-Xe(k))^2+(Yw(k)-Ye(k))^2 ) > 0.04
   Xe(k+1)= Ve*k*h;
   Ye(k+1)= H + 16*sin(2*pi*Xe(k+1)/30);
   
   Xw(k+1)= Xw(k) + Vw*h*(Xe(k)-Xw(k))/sqrt((Xe(k)-Xw(k))^2+(Ye(k)-Yw(k))^2);
   if Yw(k)<Ye(k)
     Yw(k+1)= Yw(k) + Vw*h/sqrt(1+((Xe(k)-Xw(k))/(Ye(k)-Yw(k)))^2);
   else
     Yw(k+1)= Yw(k) - Vw*h/sqrt(1+((Xe(k)-Xw(k))/(Ye(k)-Yw(k)))^2);
   end
  
   s= s+sqrt((Xw(k+1)-Xw(k))^2+(Yw(k+1)-Yw(k))^2);
   
   Wx(1)= Xw(k);
   Wx(2)= Xw(k+1);
   Wy(1)= Yw(k);
   Wy(2)= Yw(k+1);
   Ex(1)= Xe(k);
   Ex(2)= Xe(k+1);
   Ey(1)= Ye(k);
   Ey(2)= Ye(k+1);
 
   plot(Xe(k),Ye(k),'wo')
   plot(Xw(k),Yw(k),'w.')
   plot(Xe(k+1),Ye(k+1),'bo')
   plot(Xw(k+1),Yw(k+1),'r.')
   plot(Ex,Ey,'b')
   for rp=0:10
      plot(Wx,Wy,'y')
      plot(Wx,Wy,'w')
      plot(Wx,Wy,'r')
   end
 
   k= k+1; tk= tk+h;
end

plot(Xe(k),Ye(k),'ro')
plot(Xe(k),Ye(k),'y*')
text(Xe(k)-1,Ye(k)-8,'哐 !!')
sprintf(' k = %d , tk = %7.4f\n',k-1,tk)
ans= sprintf(' X=%8.4f, Y=%8.4f, T=%8.4f\n',Xe(k),Ye(k),s/Vw)
text(10,10,ans)
hold off
pause
close
clear all