exercise 5

题目2.6：Use the Euler method to calculate cannon shell trajectories ignoring both air drag and effect of air density (actually,ignoring the former automatically rules out the latter).Compare your results with those in Figure 2.4,and with the exact solution.

一、摘要
用Euler法解二元二阶常微分方程组

二、背景介绍
对于二阶的牛顿方程：

d2x/dt2 = 0，d2y/dt2 = -g

我们可以设速度Vx = dx/dt , Vy = dy/dt则可以把方程组降阶为4元一阶
常微分方程组：

dx/dt = Vx

dVx/dt = 0

dy/dt = Vy

dVy/dt = -g

然后便可用Euler法解这个方程组

三、正文

import pylab
from math import sqrt
from math import cos
from math import sin
from math import pi
class Euler2Dsimulation:
def init(self,speed,theta,timestep=0.01):
self._timestep=timestep
self._x=[0]
self._y=[0]
self._vx=[speed*cos(theta)]
self._vy=[speed*sin(theta)]
self.g=9.8
self.calculate()
self.plot()
def _accelerate(self,x,y,vx,vy):
ansax=0
ansay=-self.g
return (ansax,ansay)
def calculate(self):
while self._y[-1]>=0:
self._x.append(self._x[-1]+self._vx[-1]*self._timestep)
self._y.append(self._y[-1]+self._vy[-1]*self._timestep)
(ax,ay)=self._accelerate(self._x[-1],self._y[-1],self._vx[-1],self._vy[-1])
self._vx.append(self._vx[-1]+self._vx[-1]*self._timestep*ax)
self._vy.append(self._vy[-1]+self._timestep*ay)
def plot(self):
pylab.plot(self._x, self._y)
pylab.xlabel('x ($m$)')
pylab.ylabel('y ($m$)')
Euler2Dsimulation(700,pi/4)
Euler2Dsimulation(700,pi/6)
Euler2Dsimulation(700,pi/3)
pylab.show()

四、致谢

感谢杜威同学的帮助。感谢蔡浩老师。