第五次作业

习题2.6

---

摘要
对于习题2.6，本文的主要工作是用Euler法求解二元二阶常微分方程组.
背景知识
本题的数值模拟源于以下微分方程组

$$\frac{\mathrm{d}^2 x}{\mathrm{d} t^2}=0, \frac{\mathrm{d}^2 y}{\mathrm{d} t^2}=-g$$
令
$$v_{x}=\frac{\mathrm{d} x}{\mathrm{d} t},v_{y}=\frac{\mathrm{d} y}{\mathrm{d} t}$$
则可用以下代码进行模拟，此处选取了30,60,90三个角度进行模拟
数值模拟代码

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.001):
        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()

结果展示