第五次作业

计算物理 郭忠智2014301020087

---

一、摘要

本次作业完成了计算物理课本上第二章炮弹发射抛物线运动的近似计算，用Euler Method进行数值近似运算，并用解析的方法求得常微分方程的解析解，将两种方法所得的结果画图、以及计算出物体的落地点水平位移及最大高度，将之进行比较，可以看到，用Euler Method可以求得较好的近似结果。

二、背景

2.6 Use the Euler Method to calculate cannon shell trajectories ignoring both air drag and the 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.

三、正文

（一）题目相当于解常微分方程：

$${d^2x\over{dt^2}}=0 \tag1$$
$${d^2y\over{dt^2}}=-g\tag2$$
引入速度$v_x,v_y$有：

$$\left\{ \begin{aligned} {dx\over{dt^2}}=v_x\\ {dv_x\over{dt}}=0\\ \tag3 {dy\over{dt}}=v_y\\ {dv_y\over{dt}}=-g \end{aligned} \right.$$
应用欧拉方法，得到：

$$\left\{ \begin{aligned} x_{i+1}=x_i+v_{x,i}\Delta t\\ v_{x,i+1}=v_{x,i}\\ \tag4 y_{i+1}=y_i+v_{y,i}\Delta t\\ v_{x,i+1}=v_{x,i}-g\Delta t \end{aligned} \right.$$
用解析的求法得到（1）和（2）的解析解为：
$$x(t)=x_0+v_{x,0}t \tag4$$
$$y(t)==-{1\over2}gt^2+v_{y,0}t+y_0$$
（二）用Python编写代码，计算欧拉方法给出的结果：

import numpy
import math
import pylab as pl
class cannon_shell:
    def __init__(self, x_0, y_0, initial_velocity, theta_0, \
                 time_step, total_time):
        self.theta = theta_0
        self.x = [x_0]
        self.y = [y_0]
        self.v_x = [math.cos(theta_0 * math.pi / 180) * initial_velocity]
        self.v_y = [math.sin(theta_0 * math.pi / 180) * initial_velocity]
        self.g = 9.8
        self.dt = time_step
        self.time = total_time
    def run(self):
        for i in range(self.time):
            self.v_x.append(self.v_x[i])
            self.x.append(self.x[i] + self.dt * self.v_x[i])
            self.v_y.append(self.v_y[i] - self.g * self.dt)
            self.y.append(self.y[i] + self.dt * self.v_y[i])
            if self.y[i] >= 0 and self.y[i + 1] <= 0:
                break
        r = - self.y[-2] / self.y[-1]
        self.x[-1] = (self.x[-2] + r * self.x[-1]) / (r + 1)
        for j in range(len(self.x)):
            self.x[j] = self.x[j] / 1000
        for j in range(len(self.y)):
            self.y[j] = self.y[j] / 1000
        print "在初射角为", str(self.theta), "度时:"
        print "炮弹在水平方向的位移为：", self.x[-1], 'km'
        print "炮弹在竖直方向的最大高度为：", max(self.y), 'km'
    def show_results(self):
        font = {'family': 'serif',
                'color':  'darkred',
                'weight': 'normal',
                'size': 16,
        }
        pl.plot(self.x, self.y)
        pl.title('Trajectory of cannon shell', fontdict = font)
        pl.ylim(0, 20)
        pl.xlabel('x ($km$)')
        pl.ylabel('y ($km$)')
        pl.text(self.x[450], self.y[500],\
                str(self.theta) +'$\degree$', fontdict = font)
        pl.text(35, 18,\
                'Euler Method', fontdict = font)
a = cannon_shell(0 , 0, 700, 30, 0.1, 1500)
a.run()
a.show_results()
a = cannon_shell(0 , 0, 700, 35, 0.1, 1500)
a.run()
a.show_results()
a = cannon_shell(0 , 0, 700, 40, 0.1, 1500)
a.run()
a.show_results()
a = cannon_shell(0 , 0, 700, 45, 0.1, 1500)
a.run()
a.show_results()
a = cannon_shell(0 , 0, 700, 50, 0.1, 1500)
a.run()
a.show_results()
a = cannon_shell(0 , 0, 700, 55, 0.1, 1500)
a.run()
a.show_results()
pl.show()

Python编写代码，给出解析解的计算结果，代码如下：

#exact solution
class exact_solution:
    def __init__(self, x_0, y_0, initial_velocity, theta_0):
        self.theta = theta_0
        self.x = [x_0]
        self.y = [y_0]
        self.v_x = math.cos(theta_0 * math.pi / 180) * initial_velocity
        self.v_y = math.sin(theta_0 * math.pi / 180) * initial_velocity
        self.t = numpy.linspace(0, 150)
    def run(self):
        for i in self.t:
            self.x.append(exact_solution.xshift(self, i))
            self.y.append(exact_solution.yshift(self, i))
        r = - self.y[-2] / self.y[-1]
        self.x[-1] = (self.x[-2] + r * self.x[-1]) / (r + 1)
        for j in range(len(self.x)):
            self.x[j] = self.x[j] / 1000
        for j in range(len(self.y)):
            self.y[j] = self.y[j] / 1000
        print "在初射角为", str(self.theta), "度时:"
        print "炮弹在水平方向的位移为：", self.x[-1], 'km'
        print "炮弹在竖直方向的最大高度为：", max(self.y), 'km'
    def show_results(self):
        font = {'family': 'serif',
                'color':  'darkred',
                'weight': 'normal',
                'size': 16,
        }
        pl.plot(self.x, self.y)
        pl.title('Trajectory of cannon shell', fontdict = font)
        pl.ylim(0, 20)
        pl.xlim(0, 60)
        pl.xlabel('x ($km$)')
        pl.ylabel('y ($km$)')
        pl.text(self.x[15], self.y[15],\
                str(self.theta) +'$\degree$', fontdict = font)
        pl.text(35, 18,\
                'Exact Solution', fontdict = font)        
    def xshift(self, t):
        return self.v_x * t + self.x[0]
    def yshift(self, t):
        return -(9.8 * t ** 2) / 2 + self.v_y * t + self.y[0]
b = exact_solution(0, 0, 700, 30)
b.run()
b.show_results()
b = exact_solution(0, 0, 700, 35)
b.run()
b.show_results()  
b = exact_solution(0, 0, 700, 40)
b.run()
b.show_results() 
b = exact_solution(0, 0, 700, 45)
b.run()
b.show_results()   
b = exact_solution(0, 0, 700, 50)
b.run()
b.show_results()
b = exact_solution(0, 0, 700, 55)
b.run()
b.show_results()
pl.show()

四、结论

欧拉方法给出的结果如下：


用解析解计算，得到的结果如下：

可以得到两种方法得到的水平位移结果对比：

初射角度	欧拉方法	精确解
30	43.3618746816 km	43.2838451308 km
35	47.0419549851 km	46.972669456 km
40	49.2939993236 km	49.240204593 km
45	50.0494911508 km	49.999920421 km
50	49.2853734504 km	49.2330367646 km
55	47.024775564 km	46.9790669971 km

可以看到用两种方法得到的结果有一定误差，这是用数值算法的必然结果，可以根据实际需要，进行数据的取舍，如果对所考虑问题的精度造成的影响不大，则可以选择欧拉方法进行近似计算。