ep3_p4.m

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

%  Ep3_p4: << 实验三 >>  改进的Euler法 (预报-校正法) --- 图形演示版
%  Designed by FGH

h= 0.001;
H= 120; Vw= 450; Ve= 90;

ex= [-0.6 -0.6 0.35 0.8];
ey= [0 1 1 0];
clf
axis([-5 30 -10 130])
hold on
title('导 弹 打 敌 舰')
plot(ex,H+ey,'b')
plot(ex,H-ey,'b')
plot(0,0,'r.')
pause

tk= 0; k= 0;
% (3.28)
xk= 0; yk= 0;
while yk < H
   tk1= tk+h;
%  (3.26) 
   xsk1= xk + Vw*h*(Ve*tk-xk)/sqrt((Ve*tk-xk)^2+(H-yk)^2);
%  (3.27)   
   ysk1= yk + Vw*h/sqrt(1+((Ve*tk-xk)/(H-yk))^2);
%  (3.24)
   xk1= 0.5*(xsk1 + xk + Vw*h/sqrt(1+((H-ysk1)/(Ve*tk1-xsk1))^2));
%  (3.25)
   yk1= 0.5*(ysk1 + yk + Vw*h/sqrt(1+((Ve*tk1-xsk1)/(H-ysk1))^2));
   
   Px= [xk xk1];
   Py= [yk yk1];
   Mx= Px;
   My= [H H];
   plot(xk+ex,H+ey,'w')
   plot(xk+ex,H-ey,'w')
   plot(xk,yk,'w.')
   
%   pause
   plot(xk1+ex,H+ey,'b')
   plot(xk1+ex,H-ey,'b')
   plot(xk1,yk1,'r.')
   plot(Mx,My,'b')
   plot(Px,Py,'r')
   
   
   for rp=0:10
      plot(Px,Py,'y')
      plot(Px,Py,'w')
      plot(Px,Py,'r')
   end
   
   k= k+1;
   tk= tk+h; xk= xk1; yk= yk1;
end
plot(xk,H,'ro')
plot(xk,H,'y*')
text(xk+1.5,H-1,'哐 !!')
sprintf(' k = %d , tk = %7.4f\n',k,tk)
ans= sprintf(' L=%8.4f , T=%8.4f\n',xk,xk/Ve)
text(10,10,ans)
hold off
pause
close
clear all